library(raster)

library(plyr)
library(dplyr)

library(ggplot2)
library(viridisLite)
library(colorspace)
library(RColorBrewer)

library(sf)
library(stringr)

library(ggpubr)
library(kableExtra)

memory.limit(size = 160000)
## [1] Inf

#display.brewer.pal(n = 11, name ="RdYlGn")
colorvec <- brewer.pal(n = 11, name ="RdYlGn")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

1 Data

land_mask <- raster("Masks/land_sea_mask_1degree.nc4") 
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land


elev <- raster("Worldclim/wc2.1_10m/wc2.1_10m_elev.tif") 
elev.df <- elev %>% projectRaster(to = land_mask) %>% as.data.frame(xy = T) %>% setNames(c("lon","lat", "z"))


ipcc_regions <- shapefile("Masks/IPCC-WGI-reference-regions-v4.shp") %>% spTransform(crs("EPSG:4326")) 
ipcc_regions.raster <- ipcc_regions %>% rasterize(land_mask)
ipcc_regions.df <- as.data.frame(ipcc_regions.raster, xy = T) %>% setNames(c("lon","lat","Continent","Type","Name","Acronym"))
tas.all_annual <- read.table("CMIP6/tas/tas.all_annual.txt")
pr.all_annual <- read.table("CMIP6/pr/pr.all_annual.txt")
sfcWind.all_annual <- read.table("CMIP6/sfcWind/sfcWind.all_annual.txt")
hfls.all_annual <- read.table("CMIP6/hfls/hfls.all_annual.txt")
hfss.all_annual <- read.table("CMIP6/hfss/hfss.all_annual.txt")
hurs.all_annual <- read.table("CMIP6/hurs/hurs.all_annual.txt")

cmip6_annual <- tas.all_annual %>% merge(pr.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_annual, "cmip6_annual.txt")
tas.all_january <- read.table("CMIP6/tas/tas.all_january.txt")
pr.all_january <- read.table("CMIP6/pr/pr.all_january.txt")
sfcWind.all_january <- read.table("CMIP6/sfcWind/sfcWind.all_january.txt")
hfls.all_january <- read.table("CMIP6/hfls/hfls.all_january.txt")
hfss.all_january <- read.table("CMIP6/hfss/hfss.all_january.txt")
hurs.all_january <- read.table("CMIP6/hurs/hurs.all_january.txt")

cmip6_january <- tas.all_january %>% merge(pr.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_january, "cmip6_january.txt")
tas.all_february <- read.table("CMIP6/tas/tas.all_february.txt")
pr.all_february <- read.table("CMIP6/pr/pr.all_february.txt")
sfcWind.all_february <- read.table("CMIP6/sfcWind/sfcWind.all_february.txt")
hfls.all_february <- read.table("CMIP6/hfls/hfls.all_february.txt")
hfss.all_february <- read.table("CMIP6/hfss/hfss.all_february.txt")
hurs.all_february <- read.table("CMIP6/hurs/hurs.all_february.txt")

cmip6_february <- tas.all_february %>% merge(pr.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_february, "cmip6_february.txt")
tas.all_march <- read.table("CMIP6/tas/tas.all_march.txt")
pr.all_march <- read.table("CMIP6/pr/pr.all_march.txt")
sfcWind.all_march <- read.table("CMIP6/sfcWind/sfcWind.all_march.txt")
hfls.all_march <- read.table("CMIP6/hfls/hfls.all_march.txt")
hfss.all_march <- read.table("CMIP6/hfss/hfss.all_march.txt")
hurs.all_march <- read.table("CMIP6/hurs/hurs.all_march.txt")

cmip6_march <- tas.all_march %>% merge(pr.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_march, "cmip6_march.txt")
tas.all_april <- read.table("CMIP6/tas/tas.all_april.txt")
pr.all_april <- read.table("CMIP6/pr/pr.all_april.txt")
sfcWind.all_april <- read.table("CMIP6/sfcWind/sfcWind.all_april.txt")
hfls.all_april <- read.table("CMIP6/hfls/hfls.all_april.txt")
hfss.all_april <- read.table("CMIP6/hfss/hfss.all_april.txt")
hurs.all_april <- read.table("CMIP6/hurs/hurs.all_april.txt")

cmip6_april <- tas.all_april %>% merge(pr.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_april, "cmip6_april.txt")
tas.all_may <- read.table("CMIP6/tas/tas.all_may.txt")
pr.all_may <- read.table("CMIP6/pr/pr.all_may.txt")
sfcWind.all_may <- read.table("CMIP6/sfcWind/sfcWind.all_may.txt")
hfls.all_may <- read.table("CMIP6/hfls/hfls.all_may.txt")
hfss.all_may <- read.table("CMIP6/hfss/hfss.all_may.txt")
hurs.all_may <- read.table("CMIP6/hurs/hurs.all_may.txt")

cmip6_may <- tas.all_may %>% merge(pr.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_may, "cmip6_may.txt")
tas.all_june <- read.table("CMIP6/tas/tas.all_june.txt")
pr.all_june <- read.table("CMIP6/pr/pr.all_june.txt")
sfcWind.all_june <- read.table("CMIP6/sfcWind/sfcWind.all_june.txt")
hfls.all_june <- read.table("CMIP6/hfls/hfls.all_june.txt")
hfss.all_june <- read.table("CMIP6/hfss/hfss.all_june.txt")
hurs.all_june <- read.table("CMIP6/hurs/hurs.all_june.txt")

cmip6_june <- tas.all_june %>% merge(pr.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_june, "cmip6_june.txt")
tas.all_july <- read.table("CMIP6/tas/tas.all_july.txt")
pr.all_july <- read.table("CMIP6/pr/pr.all_july.txt")
sfcWind.all_july <- read.table("CMIP6/sfcWind/sfcWind.all_july.txt")
hfls.all_july <- read.table("CMIP6/hfls/hfls.all_july.txt")
hfss.all_july <- read.table("CMIP6/hfss/hfss.all_july.txt")
hurs.all_july <- read.table("CMIP6/hurs/hurs.all_july.txt")

cmip6_july <- tas.all_july %>% merge(pr.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_july, "cmip6_july.txt")
tas.all_august <- read.table("CMIP6/tas/tas.all_august.txt")
pr.all_august <- read.table("CMIP6/pr/pr.all_august.txt")
sfcWind.all_august <- read.table("CMIP6/sfcWind/sfcWind.all_august.txt")
hfls.all_august <- read.table("CMIP6/hfls/hfls.all_august.txt")
hfss.all_august <- read.table("CMIP6/hfss/hfss.all_august.txt")
hurs.all_august <- read.table("CMIP6/hurs/hurs.all_august.txt")

cmip6_august <- tas.all_august %>% merge(pr.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_august, "cmip6_august.txt")
tas.all_september <- read.table("CMIP6/tas/tas.all_september.txt")
pr.all_september <- read.table("CMIP6/pr/pr.all_september.txt")
sfcWind.all_september <- read.table("CMIP6/sfcWind/sfcWind.all_september.txt")
hfls.all_september <- read.table("CMIP6/hfls/hfls.all_september.txt")
hfss.all_september <- read.table("CMIP6/hfss/hfss.all_september.txt")
hurs.all_september <- read.table("CMIP6/hurs/hurs.all_september.txt")

cmip6_september <- tas.all_september %>% merge(pr.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_september, "cmip6_september.txt")
tas.all_october <- read.table("CMIP6/tas/tas.all_october.txt")
pr.all_october <- read.table("CMIP6/pr/pr.all_october.txt")
sfcWind.all_october <- read.table("CMIP6/sfcWind/sfcWind.all_october.txt")
hfls.all_october <- read.table("CMIP6/hfls/hfls.all_october.txt")
hfss.all_october <- read.table("CMIP6/hfss/hfss.all_october.txt")
hurs.all_october <- read.table("CMIP6/hurs/hurs.all_october.txt")

cmip6_october <- tas.all_october %>% merge(pr.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_october, "cmip6_october.txt")
tas.all_november <- read.table("CMIP6/tas/tas.all_november.txt")
pr.all_november <- read.table("CMIP6/pr/pr.all_november.txt")
sfcWind.all_november <- read.table("CMIP6/sfcWind/sfcWind.all_november.txt")
hfls.all_november <- read.table("CMIP6/hfls/hfls.all_november.txt")
hfss.all_november <- read.table("CMIP6/hfss/hfss.all_november.txt")
hurs.all_november <- read.table("CMIP6/hurs/hurs.all_november.txt")

cmip6_november <- tas.all_november %>% merge(pr.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_november, "cmip6_november.txt")
tas.all_december <- read.table("CMIP6/tas/tas.all_december.txt")
pr.all_december <- read.table("CMIP6/pr/pr.all_december.txt")
sfcWind.all_december <- read.table("CMIP6/sfcWind/sfcWind.all_december.txt")
hfls.all_december <- read.table("CMIP6/hfls/hfls.all_december.txt")
hfss.all_december <- read.table("CMIP6/hfss/hfss.all_december.txt")
hurs.all_december <- read.table("CMIP6/hurs/hurs.all_december.txt")

cmip6_december <- tas.all_december %>% merge(pr.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_december, "cmip6_december.txt")

2 Aridity Index

2.1 De Martonne

cmip6 <- read.table("cmip6_annual.txt")

breaks.martonne <- c(0,5,10,20,"30","40",Inf)
cat.martonne <- c("[0,5]" = "Desert", "(5,10]"= "Arid", "(10,20]" = "Semi-arid","(20,30]" = "Temperate", "(30,40]" = "Humid", "(40,Inf]" = "Forest")
col.martonne <- c("Desert" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Temperate" = colorvec[8], "Humid"= colorvec[10],"Forest"= colorvec[11])

cmip6$AIm <- with(cmip6,pr*60*60*24*365/(tas-273.15+10)) # pr in mm/y, tas in Ceslsius


cmip6$cat.AIm <- cmip6$AIm %>% cut(breaks.martonne, include.lowest = T) %>% revalue(cat.martonne)


map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1970_2000" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AIm))+
    scale_fill_manual(values = col.martonne, na.translate = F)+
    labs(title = paste(i, "1970-2000", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

2.2 UNEP

ea from Rh:

ea = 6.108(RH/100)exp(17.27*T/(T+237.3)) Refs https://www.nature.com/articles/s41598-023-40499-6 : Allen and Tetens

id_vars <- c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")

cmip6_annual <- mutate(read.table("cmip6_annual.txt"),
                       t = tas - 273.15,
                       Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                       P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                       es = 0.61078*exp(17.27*t/(t+237.3)),
                       ea = (hurs/100)*es,
                       delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                       gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                       ET0_annual = ((0.408*delta*Rn+gamma*900/((t)+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind)))) 

cmip6_01 <- mutate(read.table("cmip6_january.txt"),
                   t = tas - 273.15,
                   pr_01 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_01 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind)))) # Evapotranspiration, G considered to be negligible. in mm/day. Wind speed at 10 m is converted to wind speed at 2 m by multiplying by 0.748

cmip6_02 <- mutate(read.table("cmip6_february.txt"),
                   t = tas - 273.15,
                   pr_02 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_02 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_03 <- mutate(read.table("cmip6_march.txt"),
                   t = tas - 273.15,
                   pr_03 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_03 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_04 <- mutate(read.table("cmip6_april.txt"),
                   t = tas - 273.15,
                   pr_04 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_04 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_05 <- mutate(read.table("cmip6_may.txt"),
                   t = tas - 273.15,
                   pr_05 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_05 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_06 <- mutate(read.table("cmip6_june.txt"),
                   t = tas - 273.15,
                   pr_06 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_06 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_07 <- mutate(read.table("cmip6_july.txt"),
                   t = tas - 273.15,
                   pr_07 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_07 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_08 <- mutate(read.table("cmip6_august.txt"),
                   t = tas - 273.15,
                   pr_08 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_08 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_09 <- mutate(read.table("cmip6_september.txt"),
                   t = tas - 273.15,
                   pr_09 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_09 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_10 <- mutate(read.table("cmip6_october.txt"),
                   t = tas - 273.15,
                   pr_10 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_10 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_11 <- mutate(read.table("cmip6_november.txt"),
                   t = tas - 273.15,
                   pr_11 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_11 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_12 <- mutate(read.table("cmip6_december.txt"),
                   t = tas - 273.15,
                   pr_12 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_12 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

et0_monthly <- select(cmip6_annual, id_vars) %>%
  merge(select(cmip6_01, c(id_vars, "ET0_01", "pr_01")), by = id_vars) %>%
  merge(select(cmip6_02, c(id_vars, "ET0_02", "pr_02")), by = id_vars) %>%
  merge(select(cmip6_03, c(id_vars, "ET0_03", "pr_03")), by = id_vars) %>%
  merge(select(cmip6_04, c(id_vars, "ET0_04", "pr_04")), by = id_vars) %>%
  merge(select(cmip6_05, c(id_vars, "ET0_05", "pr_05")), by = id_vars) %>%
  merge(select(cmip6_06, c(id_vars, "ET0_06", "pr_06")), by = id_vars) %>%
  merge(select(cmip6_07, c(id_vars, "ET0_07", "pr_07")), by = id_vars) %>%
  merge(select(cmip6_08, c(id_vars, "ET0_08", "pr_08")), by = id_vars) %>%
  merge(select(cmip6_09, c(id_vars, "ET0_09", "pr_09")), by = id_vars) %>%
  merge(select(cmip6_10, c(id_vars, "ET0_10", "pr_10")), by = id_vars) %>%   
  merge(select(cmip6_11, c(id_vars, "ET0_11", "pr_11")), by = id_vars) %>%
  merge(select(cmip6_12, c(id_vars, "ET0_12", "pr_12")), by = id_vars)

et0_monthly$spr <- with(et0_monthly, pr_01*31 + pr_02 * 28 + pr_03 * 31 + pr_04 * 30 + pr_05 * 31 + pr_06 * 30 +
                          pr_07 * 31 + pr_08 * 31 + pr_09 * 30 + pr_10 * 31 + pr_11 * 30 + pr_12 * 31)
et0_monthly$sET0 <- with(et0_monthly, ET0_01*31 + ET0_02 * 28 + ET0_03 * 31 + ET0_04 * 30 + ET0_05 * 31 + ET0_06 * 30 +
                           ET0_07 * 31 + ET0_08 * 31 + ET0_09 * 30 + ET0_10 * 31 + ET0_11 * 30 + ET0_12 * 31)

cmip6 <- merge(cmip6_annual, et0_monthly, by = id_vars)
cmip6$AI <- with(cmip6, abs(spr/sET0))  
cmip6$AI_annual <- with(cmip6, abs((pr*60*60*24*365)/(ET0_annual*365)))

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

cmip6$cat.AI <- cmip6$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6$cat.AI[which(cmip6$sET0 < 400)] <- "Cold"
cmip6$cat.AI %>% unique()

cmip6$cat.AI <- cmip6$AI %>% cut(breaks = breaks.unesco) %>% 
  revalue(cat.unesco)


cmip6$cat.AI_annual <- cmip6$AI_annual %>% cut(breaks = breaks.unesco) %>% 
  revalue(cat.unesco)

cmip6$cat.AI_annual[which(cmip6$ET0_annual*365 < 400)] <- "Cold"

write.table(cmip6, "cmip6_AI.txt")

3 Maps AI by model, all % calculated without Antarctica

cmip6 <- read.table("cmip6_AI.txt")

breaks.unesco <- c(-Inf,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(-Inf,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid")

3.1 1850-1880

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1850_1880" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "1850-1880", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")


sdmm <- cmip6 %>% subset(period == "1850_1880"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))


tab <- cmip6 %>% subset(period == "1850_1880"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)


knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 5.1 9.0 13.2 6.1 33.4 38.0 28.3 0.3
CESM 5.2 10.8 14.7 5.0 35.7 35.6 28.5 0.3
CMCC 5.6 7.2 9.5 5.3 27.6 42.1 30.1 0.3
CMCC-ESM2 6.0 8.2 7.5 2.8 24.5 24.0 51.2 0.3
CNRM 5.7 12.1 11.4 6.6 35.8 38.0 25.9 0.2
EC-Earth3 8.6 9.0 11.0 4.5 33.1 36.7 30.1 0.2
FGOALS 5.6 12.8 15.9 7.2 41.5 29.9 28.4 0.2
GFDL-ESM4 5.6 9.3 8.7 4.1 27.7 37.0 34.9 0.2
INM 4.4 8.2 13.4 5.8 31.8 40.6 27.4 0.3
INM-CM5 3.5 8.8 11.5 5.6 29.4 43.2 27.1 0.3
MPI 8.5 10.9 10.1 4.4 33.9 34.5 31.4 0.2
MRI 6.3 10.6 8.7 3.8 29.4 39.6 30.9 0.3
NorESM-2-MM 5.1 10.4 18.4 5.9 39.8 35.4 24.4 0.3
mean 5.8 9.8 11.8 5.2 32.6 36.5 30.7 0.3
sd 1.4 1.6 3.2 1.2 7.4 5.1 6.7 0.1

3.2 1970-2000

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1970_2000" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "1850-1880", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

library(ggrepel) 

pie <- cmip6 %>% subset(period == "1970_2000" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  select(c("source", "cat.AI")) %>%
  reshape2::melt(id.vars = c("source")) %>%
  group_by(source, value) %>%
  summarise(count = n()) %>%
  mutate(perc = round(count/sum(count)*100, 1),
         lab.y = rev(cumsum(rev(count))) - count*0.5) %>% 
  ungroup() 

pie_list <- list()

for(i in unique(pie$source)){
  p <- ggplot(subset(pie, source == i))+
    geom_bar(aes(x = "", y = count, fill = value), width = 1, stat = "identity")+
    coord_polar("y", start = 0)+
    scale_fill_manual(values = col.cat)+
    geom_label_repel(aes(x = "", y = lab.y, label = perc), nudge_x = 0.5)+
    labs(fill = "", title = i)+
    theme_void()+theme(legend.position = "right")
  pie_list[[i]] <- p
}

ggpubr::ggarrange(plotlist = pie_list, nrow = 4, ncol = 4, common.legend = T, legend = "bottom")

pie <- cmip6 %>% subset(period == "1970_2000" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  select(c("cat.AI")) %>%
  table() %>% as.data.frame() %>%
  mutate(perc = round(Freq/sum(Freq)*100, 1),
         lab.y = rev(cumsum(rev(Freq))) - Freq*0.5) 

p <- ggplot(pie)+
    geom_bar(aes(x = "", y = Freq, fill = cat.AI), width = 1, stat = "identity")+
    coord_polar("y", start = 0)+
    scale_fill_manual(values = col.cat)+
    geom_label_repel(aes(x = "", y = lab.y, label = perc), nudge_x = 0.5)+
    labs(fill = "", title = "Multimodel CMIP6")+
    theme_void()+theme(legend.position = "right")

print(p)


sdmm <- cmip6 %>% subset(period == "1970_2000"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))


tab <- cmip6 %>% subset(period == "1970_2000"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)


knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 5.2 8.6 13.2 6.3 33.3 37.4 29.1 0.3
CESM 5.3 10.7 14.1 5.2 35.3 35.0 29.5 0.3
CMCC 5.1 7.4 10.1 5.4 28.0 42.2 29.5 0.3
CMCC-ESM2 6.4 7.9 7.7 2.5 24.5 24.4 50.8 0.3
CNRM 5.8 12.2 11.3 6.5 35.8 38.3 25.7 0.2
EC-Earth3 8.4 9.5 10.8 4.4 33.1 35.7 31.0 0.2
FGOALS 5.9 10.9 10.0 4.9 31.7 25.1 42.9 0.2
GFDL-ESM4 5.8 8.9 8.9 4.1 27.7 35.3 36.8 0.2
INM 4.4 8.4 13.3 5.7 31.8 40.1 27.8 0.3
INM-CM5 3.4 8.5 12.0 5.6 29.5 43.0 27.3 0.3
MPI 8.7 11.1 10.2 4.3 34.3 34.0 31.5 0.2
MRI 6.1 10.8 9.2 3.9 30.0 38.9 30.9 0.3
NorESM-2-MM 5.2 10.1 17.1 5.6 38.0 33.4 28.3 0.3
mean 5.8 9.6 11.4 5.0 31.8 35.6 32.4 0.3
sd 1.4 1.5 2.5 1.1 6.5 5.7 7.1 0.1

3.3 1985-2010

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1985_2015" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "1985-2010", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")


sdmm <- cmip6 %>% subset(period == "1985_2015"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))


tab <- cmip6 %>% subset(period == "1985_2015"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)


knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 5.3 8.4 13.8 6.6 34.1 37.2 28.3 0.3
CESM 4.8 10.8 14.8 5.4 35.8 35.4 28.5 0.3
CMCC 5.2 7.6 10.4 5.8 29.0 42.9 27.9 0.3
CMCC-ESM2 6.2 8.5 7.7 2.5 24.9 24.5 50.3 0.3
CNRM 5.8 12.1 11.3 6.2 35.4 39.2 25.2 0.2
EC-Earth3 8.3 9.8 10.7 4.4 33.2 36.7 29.9 0.2
FGOALS 5.8 11.4 10.4 4.9 32.5 25.3 42.0 0.2
GFDL-ESM4 5.9 9.1 9.1 3.9 28.0 36.5 35.3 0.2
INM 4.0 7.7 11.3 5.2 28.2 36.5 35.0 0.3
INM-CM5 3.5 8.3 12.0 5.4 29.2 40.6 29.9 0.3
MPI 8.7 11.1 10.3 4.4 34.5 34.1 31.0 0.2
MRI 6.2 10.7 9.1 3.8 29.8 38.4 31.6 0.3
NorESM-2-MM 4.7 9.9 15.9 5.6 36.1 30.3 33.2 0.3
mean 5.7 9.6 11.3 4.9 31.5 35.2 32.9 0.3
sd 1.5 1.5 2.3 1.1 6.4 5.5 6.7 0.1

3.4 SSP 2-4.5

3.4.1 2030-2060

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP245" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")


sdmm <- cmip6 %>% subset(period == "2030_2060" & model == "SSP245" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2030_2060" & model == "SSP245"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 5.0 9.2 14.6 6.5 35.3 38.1 26.2 0.3
CESM 4.7 12.4 15.3 5.6 38.0 35.6 26.2 0.3
CMCC 4.4 8.4 11.0 6.1 29.9 44.3 25.5 0.3
CMCC-ESM2 6.1 9.1 8.0 2.8 26.0 25.3 48.4 0.3
CNRM 6.0 12.7 12.1 6.8 37.6 39.8 22.3 0.2
EC-Earth3 7.5 9.2 11.1 4.5 32.3 36.3 31.2 0.2
FGOALS 6.9 13.3 18.5 7.8 46.5 30.3 22.9 0.2
GFDL-ESM4 5.4 10.1 9.5 4.4 29.4 36.9 33.4 0.2
INM 4.6 8.4 13.8 5.9 32.7 40.4 26.7 0.3
INM-CM5 3.5 8.6 12.2 5.4 29.7 43.1 26.9 0.3
MPI 9.0 11.9 11.1 4.5 36.5 35.1 28.2 0.2
MRI 6.2 11.2 9.3 4.1 30.8 40.0 28.9 0.3
NorESM-2-MM 5.3 11.8 16.9 5.6 39.6 33.5 26.7 0.3
mean 4.0 7.3 8.7 3.7 23.7 25.6 48.0 2.6
sd 1.0 1.2 2.2 0.9 5.3 3.6 4.6 0.0

3.4.2 2070-2100

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP245" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2070_2100" & model == "SSP245" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2070_2100" & model == "SSP245"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 4.8 10.5 15.0 7.3 37.6 38.2 24.1 0.3
CESM 6.1 11.6 15.4 5.4 38.5 36.9 24.3 0.3
CMCC 4.3 8.0 12.3 6.4 31.0 46.6 22.1 0.3
CMCC-ESM2 6.2 9.5 8.5 2.8 27.0 26.8 46.0 0.3
CNRM 6.3 12.9 12.4 6.8 38.4 41.0 20.4 0.2
EC-Earth3 7.6 9.1 11.0 4.4 32.1 37.3 30.3 0.2
FGOALS 6.9 14.0 18.1 8.1 47.1 30.3 22.3 0.2
GFDL-ESM4 5.9 9.7 9.8 4.4 29.8 38.4 31.6 0.2
INM 4.8 8.2 14.1 5.8 32.9 40.6 26.2 0.3
INM-CM5 3.7 9.5 11.2 5.9 30.3 43.5 25.9 0.3
MPI 9.0 11.9 11.1 4.5 36.5 35.1 28.2 0.2
MRI 6.5 11.0 9.6 4.0 31.1 40.8 27.9 0.3
NorESM-2-MM 5.9 11.7 17.0 5.9 40.5 34.0 25.3 0.3
mean 4.2 7.4 8.8 3.8 24.2 26.2 47.0 2.6
sd 1.0 1.3 2.1 1.0 5.4 3.7 4.5 0.0

3.5 SSP 3-7.0

3.5.1 2030-2060

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP370" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2030_2060" & model == "SSP370" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2030_2060" & model == "SSP370"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 4.9 9.6 14.3 6.8 35.6 37.4 26.8 0.3
CESM 5.2 11.4 15.6 5.4 37.6 35.7 26.5 0.3
CMCC 4.2 8.3 11.3 5.9 29.7 44.1 25.9 0.3
CMCC-ESM2 3.3 8.4 9.0 3.6 24.3 26.9 48.6 0.3
CNRM 6.2 12.9 12.9 7.1 39.1 39.6 21.0 0.2
EC-Earth3 7.9 9.7 11.6 4.8 34.0 38.5 27.3 0.2
FGOALS 6.1 11.9 11.9 5.3 35.2 27.5 37.0 0.2
GFDL-ESM4 6.3 9.4 9.9 4.1 29.7 35.9 34.1 0.2
INM 4.6 8.3 13.5 5.9 32.3 40.3 27.1 0.3
INM-CM5 3.4 8.8 12.1 5.6 29.9 43.2 26.5 0.3
MPI 8.7 12.2 11.0 4.5 36.4 36.4 26.9 0.2
MRI 6.5 10.8 9.4 3.6 30.3 40.2 29.2 0.3
NorESM-2-MM 5.4 11.0 17.9 5.7 40.0 33.2 26.5 0.3
mean 3.9 7.1 8.6 3.7 23.3 25.6 48.5 2.6
sd 1.1 1.1 1.8 0.8 4.8 3.6 4.8 0.0

3.5.2 2070-2100

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP370" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2070_2100" & model == "SSP370" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2070_2100" & model == "SSP370"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 5.5 10.5 15.4 6.9 38.3 37.3 24.2 0.3
CESM 6.0 12.2 15.6 5.7 39.5 37.4 22.9 0.3
CMCC 5.2 9.1 12.5 7.3 34.1 44.1 21.6 0.3
CMCC-ESM2 6.9 9.9 8.1 3.0 27.9 25.9 45.9 0.3
CNRM 7.3 12.9 13.2 6.7 40.1 39.3 20.3 0.2
EC-Earth3 7.6 9.8 12.4 4.7 34.5 40.3 24.9 0.2
FGOALS 6.7 13.3 13.3 4.4 37.7 28.0 34.1 0.2
GFDL-ESM4 7.5 8.5 10.9 4.7 31.6 37.0 31.1 0.2
INM 6.0 8.2 14.7 6.1 35.0 39.2 25.5 0.3
INM-CM5 3.9 10.3 14.1 6.1 34.4 40.4 24.9 0.3
MPI 8.7 12.2 11.0 4.5 36.4 36.4 26.9 0.2
MRI 7.2 11.9 9.9 4.1 33.1 39.0 27.6 0.3
NorESM-2-MM 5.9 11.8 17.5 5.9 41.1 33.7 24.9 0.3
mean 4.5 7.5 9.0 3.8 24.8 25.6 47.0 2.6
sd 0.9 1.2 1.8 0.9 4.8 3.5 4.7 0.0

3.6 SSP 5-8.5

3.6.1 2030-2060

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP585" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2030_2060" & model == "SSP585" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2030_2060" & model == "SSP585"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 4.8 9.2 15.0 7.5 36.5 37.6 25.6 0.3
CESM 5.8 11.5 15.6 5.2 38.1 36.0 25.7 0.3
CMCC 4.4 8.0 11.5 6.5 30.4 44.8 24.5 0.3
CMCC-ESM2 6.5 9.3 8.0 2.5 26.3 25.4 48.0 0.3
CNRM 6.2 13.3 12.4 6.9 38.8 40.0 21.0 0.2
EC-Earth3 7.8 9.4 11.8 4.6 33.6 38.8 27.4 0.2
FGOALS 6.8 13.5 18.1 7.7 46.1 30.8 22.8 0.2
GFDL-ESM4 6.1 9.4 9.7 4.7 29.9 36.4 33.4 0.2
INM 4.6 8.4 14.0 6.0 33.0 40.6 26.2 0.3
INM-CM5 3.5 8.6 12.4 5.8 30.3 43.6 25.8 0.3
MPI 8.6 12.4 11.1 4.5 36.6 37.0 26.2 0.2
MRI 6.0 11.1 9.6 4.0 30.7 40.8 28.1 0.3
NorESM-2-MM 5.5 11.4 17.6 5.9 40.4 33.2 26.2 0.3
mean 4.1 7.3 8.9 3.8 24.1 26.0 47.3 2.6
sd 1.0 1.3 2.1 1.0 5.4 3.7 4.7 0.0

3.6.2 2070-2100

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP585" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2070_2100" & model == "SSP585" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2070_2100" & model == "SSP585"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(13,14), bold = T) 
source Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
CAS-ESM2 6.1 10.8 15.9 6.5 39.3 37.5 22.9 0.3
CMCC 4.4 9.3 13.5 6.7 33.9 46.7 19.2 0.3
CMCC-ESM2 5.6 10.5 8.9 2.9 27.9 28.1 43.8 0.3
CNRM 6.0 13.6 12.8 6.4 38.8 40.9 20.1 0.2
EC-Earth3 7.7 9.6 13.0 4.9 35.2 40.5 23.9 0.2
FGOALS 15.6 15.0 18.3 6.8 55.7 24.4 19.8 0.2
GFDL-ESM4 6.5 10.0 10.8 4.8 32.1 37.6 30.1 0.2
INM 4.4 8.7 14.8 6.0 33.9 41.7 24.1 0.3
INM-CM5 3.9 9.8 12.0 6.0 31.7 44.1 23.9 0.3
MPI 8.7 12.3 11.1 4.7 36.8 36.8 26.2 0.2
MRI 7.1 11.5 10.6 4.1 33.3 40.1 26.4 0.3
NorESM-2-MM 6.4 12.8 17.4 6.5 43.1 34.3 22.3 0.3
mean 4.8 7.8 9.2 3.8 25.6 26.2 45.5 2.6
sd 2.1 1.3 2.0 0.8 6.2 4.4 4.6 0.0

4 Agreement between models

4.1 Summary dataframe

4.1.1 Difference with AI ref for all SSP

test <- cmip6ext %>% subset(model %in% c("historical","SSP585") & period %in% c("1970_2000", "2070_2100") & lon == 25.5 & lat == 10.5) %>% select(c("source","lon","lat","model","period","AI","cat.AI","AI.ref"))
                                                                               
                                                                               
test2 <- cmip6s %>% subset(model == "SSP585" & period == "2070_2100" & lon == 25.5 & lat == 10.5) %>% 
  select(c("model","period","AI.mean","cat.AI","AI.ref.mean","diff.AI.mean", "nb.models", "diff.AI.mean"))

test2
cmip6 <- read.table("cmip6_AI.txt")

cmip6.ref <- cmip6 %>% filter(period == "1970_2000") %>% select(c("lon","lat","model","period","source","AI","spr","t")) %>% setNames(c("lon","lat","model","period","source","AI.ref","spr.ref","t.ref"))

cmip6ext <- merge(cmip6, select(cmip6.ref, c("lon","lat","source","AI.ref","spr.ref","t.ref")) , by = c("lon","lat","source"), all = T) %>%  mutate(diff.AI = AI - AI.ref, diff.pr = spr - spr.ref, diff.t = t - t.ref)

write.table(cmip6ext, "cmip6ext.txt")
cmip6 <- read.table("cmip6ext.txt")
cmip6s <- cmip6 %>% 
  group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, period, z) %>% 
  dplyr::summarise(t.mean = mean(t, na.rm = T), t.sd = sd(t, na.rm = T), # t in °C
                   pr.mean = mean(spr, na.rm = T) , pr.sd = sd(spr, na.rm = T),  # pr in mm/y
                   sfcWind.mean = mean(sfcWind, na.rm = T), sfcWind.sd = sd(sfcWind, na.rm = T),
                   Rn.mean = mean(Rn, na.rm = T), Rn.sd = sd(Rn, na.rm = T), # in MJ/m2/day
                   ea.mean = mean(ea, na.rm = T), ea.sd = sd(ea, na.rm = T), # in kPa
                   ET0 = mean(sET0, na.rm = T), ET0.sd = sd(sET0, na.rm = T), # in mm/y
                   AI.mean = mean(AI, na.rm = T), AI.sd = sd(AI, na.rm = T),
                   AI.ref.mean = mean(AI.ref, na.rm = T), AI.ref.sd = sd(AI.ref, na.rm = T),
                   diff.AI.mean = mean(diff.AI, na.rm = T), diff.AI.sd = sd(diff.AI, na.rm = T),
                   spr.ref.mean = mean(spr.ref, na.rm = T), spr.ref.sd = sd(spr.ref, na.rm = T),
                   diff.pr.mean = mean(diff.pr, na.rm = T), diff.pr.sd = sd(diff.pr, na.rm = T),
                   t.ref.mean = mean(t.ref, na.rm = T), t.ref.sd = sd(t.ref, na.rm =T),
                   diff.t.mean = mean(diff.t, na.rm = T), diff.t.sd = sd(diff.t, na.rm = T),
                   Hyperarid = sum(cat.AI == "Hyperarid"),
                   Arid = sum(cat.AI == "Arid"),
                   'Semi-arid' = sum(cat.AI == "Semi-arid"),
                   'Dry subhumid' = sum(cat.AI == "Dry subhumid") ,
                   'Humid' = sum(cat.AI == "Humid"),
                   'Cold' = sum(cat.AI == "Cold"),
                   'NA' = sum(is.na(cat.AI)),
                   nb.models = n(),
                    plus.AI = table(diff.AI > 0)["TRUE"],
                    minus.AI = table(diff.AI < 0)["TRUE"]) %>%
ungroup() %>% as.data.frame()

cmip6s$cat.maj <- colnames(select(cmip6s, c("Hyperarid","Arid",'Semi-arid',"Dry subhumid",'Humid',"Cold","NA")))[max.col(select(cmip6s, c("Hyperarid","Arid",'Semi-arid',"Dry subhumid",'Humid',"Cold","NA")),ties.method = "first")]
cmip6s$nb.maj <- cmip6s %>% select(c("Hyperarid","Arid",'Semi-arid',"Dry subhumid",'Humid',"Cold","NA")) %>% apply(1, max)
cmip6s$prop.maj <- with(cmip6s, nb.maj / nb.models)



breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")

cmip6s$cat.AI <- cmip6s$AI.mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6s$cat.AI[which(cmip6s$ET0 < 400)] <- "Cold"

cmip6s$cat.AI.ref <- cmip6s$AI.ref.mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6s$cat.AI.ref[which(cmip6s$ET0 < 400)] <- "Cold"

cmip6s$same.catAI <- with(cmip6s, ifelse(cat.maj == cat.AI, "same", "diff"))

cmip6s$cat.dryer <- with(cmip6s, 
          ifelse(cat.AI.ref == "Arid" & cat.AI %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
          ifelse(cat.AI.ref == "Arid" & cat.AI == "Hyperarid", "dryer",
          ifelse(cat.AI.ref == "Semi-arid" & cat.AI %in% c("Cold","Humid","Dry subhumid"), "wetter",
          ifelse(cat.AI.ref == "Semi-arid" & cat.AI %in% c("Hyperarid","Arid"), "dryer",
          ifelse(cat.AI.ref == "Dry subhumid" & cat.AI %in% c("Cold","Humid"), "wetter",
          ifelse(cat.AI.ref == "Dry subhumid" & cat.AI %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.ref == "Humid" & cat.AI == "Cold","wetter",
          ifelse(cat.AI.ref == "Humid" & cat.AI %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.ref == "Cold" & cat.AI %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.ref == "Hyperarid" & cat.AI %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
          ifelse(cat.AI.ref == cat.AI, "same",
          NA))))))))))))

cmip6sref <- cmip6s %>% filter(period == "1970_2000") %>% select(c("lon","lat","model","period","cat.maj")) %>% setNames(c("lon","lat","model","period","cat.maj.ref"))

cmip6s <- merge(cmip6s, select(cmip6sref, c("lon","lat", "model","period","cat.maj.ref")), 
                 by = c("lon","lat","model","period"), all = T) 

cmip6s$cat.maj.dryer <- with(cmip6s, 
        ifelse(cat.maj.ref == "Arid" & cat.maj %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
        ifelse(cat.maj.ref == "Arid" & cat.maj == "Hyperarid", "dryer",
        ifelse(cat.maj.ref == "Semi-arid" & cat.maj %in% c("Cold","Humid","Dry subhumid"), "wetter",
        ifelse(cat.maj.ref == "Semi-arid" & cat.maj %in% c("Hyperarid","Arid"), "dryer",
        ifelse(cat.maj.ref == "Dry subhumid" & cat.maj %in% c("Cold","Humid"), "wetter",
        ifelse(cat.maj.ref == "Dry subhumid" & cat.maj %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
        ifelse(cat.maj.ref == "Humid" & cat.maj == "Cold","wetter",
        ifelse(cat.maj.ref == "Humid" & cat.maj %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
        ifelse(cat.maj.ref == "Cold" & cat.maj %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
        ifelse(cat.maj.ref == "Hyperarid" & cat.maj %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
        ifelse(cat.maj.ref == cat.maj, "same",
        NA))))))))))))


write.table(cmip6s, "cmip6s.txt")
cmip6s <- read.table("cmip6s.txt")

cmip6s$Continent[which(cmip6s$Continent == "Europe-Africa")] <- "Mediterranea"

4.2 Compare cat mean and cat maj

cat mean = category of aridity obtained by using the multimodel mean AI cat maj = category of aridity obtained by using the category majoritarily obtained in the 13 models

4.2.1 Where are cat mean and cat maj different


cmip6s %>% group_by(model, period) %>% summarise(maj.equal.mean = table(same.catAI)[2], maj.diff.mean = table(same.catAI)[1], prop.same = round(maj.equal.mean/n()*100,1)) %>% knitr::kable(col.names = c("Model","Period","Number of gridcells in which cat.maj = cat.mean", "Number of gridcells in which cat.maj /= cat.mean", "Proportion of cells with the same category (%)")) %>% kableExtra::kable_styling(bootstrap_options = "bordered")
Model Period Number of gridcells in which cat.maj = cat.mean Number of gridcells in which cat.maj /= cat.mean Proportion of cells with the same category (%)
SSP245 2030_2060 19884 1798 89.3
SSP245 2070_2100 19890 1792 89.3
SSP370 2030_2060 19988 1694 89.8
SSP370 2070_2100 19911 1771 89.4
SSP585 2030_2060 19799 1883 88.9
SSP585 2070_2100 19796 1886 88.9
historical 1850_1880 19907 1775 89.4
historical 1970_2000 19880 1802 89.3
historical 1985_2015 19936 1746 89.6

cmip6s$cat.AI <- factor(cmip6s$cat.AI, levels =  c("Hyperarid","Arid","Semi-arid","Dry subhumid", "Humid","Cold"))
cmip6s$cat.maj <- factor(cmip6s$cat.maj, levels =  c("Hyperarid","Arid","Semi-arid","Dry subhumid", "Humid","Cold"))

ggplot(subset(cmip6s, cat.AI != "Cold"))+
  geom_col(aes(x= period, y = cat.AI, fill = cat.AI), just = -0.2, width = 0.3)+
  geom_col(aes(x = period, y = cat.maj, fill = cat.maj), just = 1.2, width = 0.3)+
  scale_fill_manual(values = col.cat, na.translate = F)+
  labs(fill = "", y ="")+
  facet_grid(rows = vars(model), switch = "y")+
  theme_minimal()+theme(axis.text.y = element_blank())




tab <- rbind(table(cmip6s$cat.AI)/length(cmip6s$cat.AI)*100, table(cmip6s$cat.maj)/length(cmip6s$cat.maj)*100) %>% t() %>% as.data.frame() %>% setNames(c("Cat mean", "Cat maj"))
knitr::kable(tab, digits =2, caption = "Proportion of gridcells in each aridity category") %>% kable_styling(bootstrap_options = "bordered")
Proportion of gridcells in each aridity category
Cat mean Cat maj
Hyperarid 3.51 4.18
Arid 6.25 7.43
Semi-arid 7.63 9.19
Dry subhumid 3.39 0.84
Humid 28.75 28.53
Cold 47.90 47.23
map_list <- list()

for(i in unique(cmip6s$model)){
  df <- subset(cmip6s, model == i)
  for(j in unique(df$period)){
    index <- paste(i, j, sep = "")
    g <- ggplot()+geom_raster(data = subset(df, period == j & same.catAI == "diff"), aes(x = lon, y = lat, fill = cat.AI))+
      borders(colour = "grey60")+
      scale_fill_manual(values =  col.cat, na.translate = F)+
      labs(title = paste(i, j , sep = " "), fill = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }
}

col1 <- ggarrange(plotlist = map_list, nrow = 9, ncol = 1, common.legend = T, legend = "bottom") %>% annotate_figure(top = "cat.mean for gridcells in which\ncat.maj and cat.mean are different")

map_list2 <- list()

for(i in unique(cmip6s$model)){
  df <- subset(cmip6s, model == i)
  for(j in unique(df$period)){
    index <- paste(i, j, sep = "")
    g <- ggplot()+geom_raster(data = subset(df, period == j & same.catAI == "diff"), aes(x = lon, y = lat, fill = cat.maj))+
      borders(colour = "grey60")+
      scale_fill_manual(values =  col.cat, na.translate = F)+
      labs(title = paste(i, j , sep = " "), fill = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list2[[index]] <- g
  }
}

col2 <- ggarrange(plotlist = map_list2, nrow = 9, ncol = 1, common.legend = T, legend = "bottom") %>% annotate_figure(top = "cat.maj for gridcells in which\ncat.maj and cat.mean are different")

ggarrange(col1,col2,ncol = 2, common.legend = T)

4.2.2 Plot cat maj and cat mean for reference period

ggpubr::ggarrange(plotlist = list(
  
  ggplot() + 
    geom_raster(data = subset(cmip6s, period == "1970_2000"), aes(x=lon, y = lat,  fill = cat.maj))+
    geom_point(data = subset(cmip6s, period == "1970_2000" & nb.maj > 10), aes(x =lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
    scale_fill_manual(values = col.cat)+
    labs(title = "Category chosen in majority by models", fill = "")+
    ylim(-55,90)+
    theme_void(base_size = 15)+
    theme(legend.position = "none"),
  
  ggplot() + 
    geom_raster(data = subset(cmip6s, period == "1970_2000"), aes(x=lon, y = lat,  fill = cat.AI))+
   geom_point(data = subset(cmip6s, period == "1970_2000" & nb.maj > 10), aes(x =lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
    scale_fill_manual(values = col.cat)+
    labs(title = "Multimodel average categories of AI", fill = "")+
    ylim(-55,90)+
    theme_void(base_size = 15)+
    theme(legend.position = "none")),
  
  ncol = 2)
map_list <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = cat.maj))+
    geom_point(data = subset(cmip6s, period == i & model == "historical" & nb.maj > 10), aes(x= lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list[[i]] <- g
}


map_list2 <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI))+
   # geom_point(data = subset(cmip6s, period == i & model == "historical" & nb.maj > 10), aes(x= lon, y = lat), size = 0.5, col = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list2[[i]] <- g
}


ggarrange(
  annotate_figure(ggarrange(plotlist = map_list, nrow = 3, common.legend = T, legend = "bottom"), top = "Majority category among models"),
  annotate_figure(ggarrange(plotlist = map_list2, nrow = 3, common.legend = T, legend = "bottom"), top = "Multimodel average category"),
  ncol = 2, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = cat.maj))+
         geom_point(data = subset(cmip6s, period == j & model == i & nb.maj > 10), aes(x= lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
      scale_fill_manual(values = col.cat, na.translate = F)+
      labs(fill = "Cat AI maj", title = paste(i, j, sep = ", "))+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}


map_list2 <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = cat.AI))+
    #     geom_point(data = subset(cmip6s, period == j & model == i & nb.maj > 10), aes(x= lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
      scale_fill_manual(values = col.cat, na.translate = F)+
      labs(fill = "Cat AI mean", title = paste(i, j, sep = ", "))+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list2[[index]] <- g
  }}

ggarrange(plotlist = map_list, nrow = 3, ncol = 2, common.legend = T, legend = "bottom")

ggarrange(plotlist = map_list2, nrow = 3, ncol = 2, common.legend = T, legend = "bottom")

4.3 Multimodel maps of cat mean


cmip6s$cat.AI <- factor(cmip6s$cat.AI, levels = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))

map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + 
    geom_tile(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI))+
    geom_point(data = subset(cmip6s,nb.maj > 10), aes(x = lon, y = lat), shape = ".", col = "grey20", size = 0.1)+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = paste("Category AI",i , sep = " "), fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + 
      geom_tile(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = cat.AI))+
      geom_point(data = subset(cmip6s,nb.maj > 10), aes(x = lon, y = lat), shape = ".", col = "grey20", size = 0.1)+
      borders(colour = "grey60")+
      scale_fill_manual(values =  col.cat, na.translate = F)+
      labs(title = paste("Category AI",index, sep = " "), fill = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

5 Maps multimodel AI

quantile(cmip6s$AI.mean, probs = seq(0,1,0.1), na.rm = T) 
##           0%          10%          20%          30%          40%          50% 
## 2.902197e-03 1.993698e-01 5.922059e-01 1.019621e+00 1.457766e+00 2.127218e+00 
##          60%          70%          80%          90%         100% 
## 3.578980e+00 1.168470e+01 2.438867e+01 4.856852e+01 5.329875e+06

ai.breaks <- c(0,0.03,0.2,0.5,0.65,1,10,20,30)
colscale <- colorvec[-c(5,6)]
map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = AI.mean))+
    borders(colour = "grey60")+
    binned_scale(aesthetics = "fill", breaks = ai.breaks, palette = function(x) colscale,
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "Aridity index")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = AI.mean))+
      borders(colour = "grey60")+
      binned_scale(aesthetics = "fill", breaks = ai.breaks, palette = function(x) colscale,
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Aridity index", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

6 AI changes over time and scenarios

v<-NULL
v2<-NULL
for (N in 5:30) {
  v<-c(v,sum(sapply(floor(2*N/3):N,function(i) choose(N,i)*(1/2)^N*2)))# proba d'avoir 2/3 de même signe à partir de N modèles, lorsque N change
  v2<-c(v2,sum(sapply(floor(2*N/3):N,function(i) choose(N,i)*(1/2)^N)))# proba d'avoir 2/3 de signe >0 à partir de N modèles, lorsque N change
}


plot(5:30,v,type="l",ylab="Proba(>2/3 models...)",xlab="Nb of models")
lines(5:30,v2,col="red")

6.1 Changes AI for summary CMIP6

6.2 Boxplot of changes AI

ggpubr::ggarrange(plotlist = list(
  ggplot(subset(cmip6s, diff.AI.mean < 1 & diff.AI.mean > -1 & model == "historical"))+
    geom_boxplot(aes(x=period, y = - diff.AI.mean, col = period), outliers = F)+
    scale_color_manual(values = c("#1A9850","#66BD63","#A6D96A"))+
    labs(x = "", y = "Difference of AI between each period and\nthe reference period 1970-2000")+
    facet_grid(cols = vars(model))+
    theme_minimal() + theme(legend.position = "none", axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)),
  ggplot(subset(cmip6s, diff.AI.mean < 1 & diff.AI.mean > -1 & model != "historical"))+
    geom_boxplot(aes(x=period, y = diff.AI.mean, col = period), outliers = F)+
    scale_color_manual(values = c("#FDAE61","#F46D43"))+
    scale_y_continuous(position = "right")+
    labs(y = "", x = "")+
    facet_grid(cols = vars(model))+
    theme_minimal() + theme(legend.position = "none", axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))),
  widths = c(1.5,3),
  ncol = 2, nrow = 1
)

6.3 Strictly wetter or dryer (reference period: 1970-2000)

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
  index <- paste(i, j, sep = ", ")
  g <- ggplot()+
  geom_point(data = subset(cmip6s, model == i & period == j & plus.AI > 10), aes(x = lon, y = lat, col = "plus"), shape = "+")+
  geom_point(data = subset(cmip6s, model == i & period == j & minus.AI > 10), aes(x = lon, y = lat, col = "minus"), shape = "-")+
    scale_color_manual(values = colorvec[c(1,11)])+
    labs(color = "", x = "", y = "", title = index)+
  theme_minimal()
  map_list[[index]] <- g
  }
}

ggarrange(plotlist = map_list, ncol = 2, nrow = 3, common.legend = T)

map_list <- list()

for(i in c("1850_1880","1985_2015")){
    index <- paste(i, j, sep = " , ")
  g <- ggplot() + geom_tile(data = subset(cmip6s, period == i & model == "historical" & diff.AI.mean > 0), aes(x=lon, y = lat,  fill = "wetter"))+
    geom_tile(data = subset(cmip6s, period == i & model == "historical" & diff.AI.mean < 0), aes(x=lon, y = lat,  fill = "dryer"))+
    geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = "+", size = 0.5)+
     geom_point(data = subset(cmip6s, period == i & model == "historical" & minus.AI > 10), aes(x=lon, y = lat), shape = "-", size = 0.5)+
    borders(colour = "grey60")+
    scale_fill_manual(values =  c("#F46D43", "#74ADD1"), na.translate = F)+
    labs(title = paste("AI changes between",i , "and 1970-2000", sep = " "), fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list[[i]] <- g
  }

ggarrange(plotlist = map_list, ncol = 2)


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & diff.AI.mean > 0), aes(x=lon, y = lat,  fill = "wetter"))+
      geom_tile(data = subset(cmip6s, period == j & model == i & diff.AI.mean < 0), aes(x=lon, y = lat,  fill = "dryer"))+
       geom_point(data = subset(cmip6s, period == j & model == i & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s, period == j & model == i & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
      borders(colour = "grey60")+
      scale_fill_manual(values = c("#F46D43", "#74ADD1"), na.translate = F)+
      labs(fill = "Compared to 1970-2000", title = i)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}



ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

6.4 Wetter or dryer, binned (reference period: 1970-2000)

ano.quantiles <- cmip6s %>% group_by(model, period) %>% summarise(q10 = quantile(diff.AI.mean, probs = 0.1, na.rm =T), 
                                                                  q25 = quantile(diff.AI.mean, probs = 0.25, na.rm =T), 
                                                                  q50 = quantile(diff.AI.mean, probs = 0.5, na.rm =T), 
                                                                  q75 = quantile(diff.AI.mean, probs = 0.75, na.rm =T), 
                                                                  q90 = quantile(diff.AI.mean, probs = 0.9, na.rm =T)) 

ano.quantiles <- cmip6s %>% summarise(q10 = quantile(diff.AI.mean, probs = 0.1, na.rm =T), 
                                      q25 = quantile(diff.AI.mean, probs = 0.25, na.rm =T), 
                                      q50 = quantile(diff.AI.mean, probs = 0.5, na.rm =T), 
                                      q75 = quantile(diff.AI.mean, probs = 0.75, na.rm =T), 
                                      q90 = quantile(diff.AI.mean, probs = 0.9, na.rm =T)) 
qbreaks <- c(-5, -1,-0.1, -0.01, 0, 0.01, 0.1, 1,5)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")

“diff.AI.mean” is the multimodel AI average minus the reference AI (period 1970-2000). More blue: AI mean is superior to AI ref, hence wetter than the ref. More red: the anomalie is negative, hence AI mean is inferior to AI ref, hence dryer.

AS OF NOW THE SAME QUANTILES ARE USED FOR ALL MAPS.

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  #quant = subset(ano.quantiles, model == "historical" & period == i)
  #b = quant[1,c(3:7)] %>% as.numeric
  b = qbreaks
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = diff.AI.mean))+
    borders(colour = "grey60")+
           geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s,  period == i & model == "historical"  & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
    binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale),
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    b = qbreaks
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = diff.AI.mean))+
      borders(colour = "grey60")+
                 geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s,  period == i & model == "historical"  & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
      binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale),
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Compared to 1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}



ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

6.5 Wetter or dryer, binned, % (reference period: 1970-2000)

cmip6s$diff.AI.mean.percent <- with(cmip6s, (AI.mean-AI.ref.mean)/AI.ref.mean*100)


pbreaks <- c(-40, -30,-20, -10, 0, 10, 20, 30, 40)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef","white","white", "#fddbc7", "#f4a582", "#d6604d", "#b2182b")

“diff.AI.mean” is the multimodel AI average minus the reference AI (period 1970-2000). More blue: AI mean is superior to AI ref, hence wetter than the ref. More red: the anomalie is negative, hence AI mean is inferior to AI ref, hence dryer.

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  #quant = subset(ano.quantiles, model == "historical" & period == i)
  #b = quant[1,c(3:7)] %>% as.numeric
  b = pbreaks
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = diff.AI.mean.percent))+
    borders(colour = "grey60")+
           geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s,  period == i & model == "historical"  & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
    binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale))+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom") %>% annotate_figure(top = "% change of AI index compared to 1970-2000")


#,guide = guide_legend(label.theme = element_text(angle = 0))

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    b = pbreaks
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = diff.AI.mean.percent))+
      borders(colour = "grey60", size = 0.5)+
                 geom_point(data = subset(cmip6s, period == j & model == i & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s,  period == j & model == i & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
      binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale))+
      labs(fill = "% change of AI index\ncompared to 1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}


ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

7 Difference between future and reference category of AI, using cat maj

7.1 Detailed changes and colors cat maj

cmip6s$change.catmaj <- with(cmip6s, ifelse(cat.maj == cat.maj.ref, 0, 1))

cmip6s$change.cat.maj <- with(cmip6s, paste(cat.maj.ref, cat.maj, sep = " to "))

cmip6s$change.cat.maj <- with(cmip6s, ifelse(cat.maj.ref == "Arid" & cat.maj %in% c("Cold","Humid"),"Arid to cold or humid",
         ifelse(cat.maj.ref == "Arid" & cat.maj %in% c("Semi-Arid","Dry subhumid"), "Arid to semi-arid or dry subhumid",
         ifelse(cat.maj.ref == "Arid" & cat.maj == "Hyperarid", "Arid to hyperarid",
         ifelse(cat.maj.ref == "Semi-arid" & cat.maj == "Dry subhumid", "Semi-arid to dry subhumid",
         ifelse(cat.maj.ref == "Semi-arid" & cat.maj == "Arid", "Semi-arid to arid",
         ifelse(cat.maj.ref == "Semi-arid" & cat.maj == "Hyperarid", "Semi-arid to hyperarid",
         ifelse(cat.maj.ref == "Semi-arid" & cat.maj %in% c("Humid","Cold"), "Semi-arid to humid or cold",
         ifelse(cat.maj.ref == "Dry subhumid" & cat.maj %in% c("Semi-arid","Arid","Hyperarid"), "Dry subhumid to arid, semi-arid or hyperarid",
         ifelse(cat.maj.ref == "Dry subhumid" & cat.maj %in% c("Humid","Cold"), "Dry subhumid to humid or cold",
         ifelse(cat.maj.ref == "Humid" & cat.maj == "Cold","Humid to cold",
         ifelse(cat.maj.ref == "Humid" & cat.maj %in% c("Dry subhumid", "Semi-arid"), "Humid to dry-subhumid or semi-arid",
         ifelse(cat.maj.ref == "Humid" & cat.maj %in% c("Arid", "Hyperarid"), "Humid to arid or hyperarid",
         ifelse(cat.maj.ref == "Cold" & cat.maj == "Humid", "Cold to humid", 
         ifelse(cat.maj.ref == "Cold" & cat.maj %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "Cold to dryland",
         ifelse(cat.maj.ref == "Hyperarid" & cat.maj == "Arid", "Hyperarid to arid",
         ifelse(cat.maj.ref == "Hyperarid" & cat.maj %in% c("Semi-arid","Dry-subhumid"),"Hyperarid to semi arid or dry subhumid",
         ifelse(cat.maj.ref == "Hyperarid" & cat.maj %in% c("Humid", "Cold"), "Hyperarid to humid or cold", NA))))))))))))))))))

cmip6s$change.cat.maj <- factor(cmip6s$change.cat.maj, 
                                levels = c("Arid to hyperarid","Semi-arid to arid", "Dry subhumid to arid, semi-arid or hyperarid", "Humid to dry-subhumid or semi-arid", "Cold to humid","Hyperarid to arid", "Arid to semi-arid or dry subhumid", "Dry subhumid to humid or cold", "Semi-arid to dry subhumid", "Semi-arid to humid or cold"))
table(cmip6s$change.cat.maj)
## 
##                            Arid to hyperarid 
##                                            0 
##                            Semi-arid to arid 
##                                            0 
## Dry subhumid to arid, semi-arid or hyperarid 
##                                            0 
##           Humid to dry-subhumid or semi-arid 
##                                            0 
##                                Cold to humid 
##                                            0 
##                            Hyperarid to arid 
##                                            0 
##            Arid to semi-arid or dry subhumid 
##                                            0 
##                Dry subhumid to humid or cold 
##                                            0 
##                    Semi-arid to dry subhumid 
##                                            0 
##                   Semi-arid to humid or cold 
##                                            0

colors_changes <- c("Hyperarid to arid" = "#E7D7B6", 
                    "Arid to semi-arid or dry subhumid" = colorvec[7],
                    "Dry subhumid to humid or cold" = colorvec[8],
                    "Semi-arid to dry subhumid" = colorvec[9],
                    "Semi-arid to humid or cold" = colorvec[10],
                    "Humid to cold" = colorvec[11],
                    "Cold to humid" = "#E7D7B6", 
                    "Humid to dry-subhumid or semi-arid" = colorvec[5],
                    "Dry subhumid to arid, semi-arid or hyperarid" = colorvec[4],
                    "Semi-arid to arid" = colorvec[3],
                    "Arid to hyperarid" = colorvec[1])

7.2 Changes towards dryer categories cat maj

7.2.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + 
    geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & cat.maj.dryer == "dryer"), 
                aes(x=lon, y = lat,  fill = change.cat.maj))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}


ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

7.2.1.1 Future


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmaj == 1 & cat.maj.dryer == "dryer"), 
                              aes(x=lon, y = lat,  fill = change.cat.maj))+
      #   geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & nb.maj > 10), aes(x= lon, y = lat, fill = "'"))+
      borders(colour = "grey60")+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      #  scale_fill_viridis_d()+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

7.3 Changes towards wetter categories cat maj

7.3.0.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & cat.maj.dryer == "wetter"), 
                              aes(x=lon, y = lat,  fill = change.cat.maj))+
    borders(colour = "grey60", ylim = c(-40,40))+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+
    ylim(-50,50)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

7.3.0.2 Future


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmaj == 1 & cat.maj.dryer == "wetter"), 
                              aes(x=lon, y = lat,  fill = change.cat.maj))+
      borders(colour = "grey60", ylim = c(-40,40))+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+ylim(-50,50)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list, nrow = 3, ncol = 2, common.legend = T, legend = "bottom")

8 Difference between future and reference category of AI, using cat mean

8.1 Detailed changes and colors

cmip6s$change.catmean <- with(cmip6s, ifelse(cat.AI == cat.AI.ref, 0, 1))

cmip6s$change.cat.mean <- with(cmip6s, paste(cat.AI.ref, cat.AI, sep = " to "))
unique(cmip6s$change.cat)
## NULL

# cmip6s$change.cat.mean <- with(cmip6s, 
#       ifelse(cat.AI.ref == "Arid" & cat.AI %in% c("Cold","Humid"),"Arid to cold or humid",
#       ifelse(cat.AI.ref == "Arid" & cat.AI %in% c("Semi-Arid","Dry subhumid"), "Arid to semi-arid or dry subhumid",
#       ifelse(cat.AI.ref == "Arid" & cat.AI == "Hyperarid", "Arid to hyperarid",
#       ifelse(cat.AI.ref == "Semi-arid" & cat.AI == "Dry subhumid", "Semi-arid to dry subhumid",
#       ifelse(cat.AI.ref == "Semi-arid" & cat.AI == "Arid", "Semi-arid to arid",
#       ifelse(cat.AI.ref == "Semi-arid" & cat.AI == "Hyperarid", "Semi-arid to hyperarid",
#       ifelse(cat.AI.ref == "Semi-arid" & cat.AI %in% c("Humid","Cold"), "Semi-arid to humid or cold",
#       ifelse(cat.AI.ref == "Dry subhumid" & cat.AI %in% c("Semi-arid","Arid","Hyperarid"), "Dry subhumid to arid, semi-arid or hyperarid",
#       ifelse(cat.AI.ref == "Dry subhumid" & cat.AI %in% c("Humid","Cold"), "Dry subhumid to humid or cold",
#       ifelse(cat.AI.ref == "Humid" & cat.AI == "Cold","Humid to cold",
#       ifelse(cat.AI.ref == "Humid" & cat.AI %in% c("Dry subhumid","Semi-arid"), "Humid to dry-subhumid or semi-arid",
#       ifelse(cat.AI.ref == "Humid" & cat.AI %in% c("Arid", "Hyperarid"), "Humid to arid or hyperarid",
#       ifelse(cat.AI.ref == "Cold" & cat.AI == "Humid", "Cold to humid", 
#       ifelse(cat.AI.ref == "Cold" & cat.AI %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "Cold to dryland",
#       ifelse(cat.AI.ref == "Hyperarid" & cat.AI == "Arid", "Hyperarid to arid",
#       ifelse(cat.AI.ref == "Hyperarid" & cat.AI %in% c("Semi-arid","Dry-subhumid"),"Hyperarid to semi arid or dry subhumid",
#       ifelse(cat.AI.ref == "Hyperarid" & cat.AI %in% c("Humid", "Cold"), "Hyperarid to humid or cold", 
#       NA))))))))))))))))))

table(cmip6s$change.cat.mean)
## 
##                 Arid to Arid            Arid to Hyperarid 
##                        11080                          639 
##            Arid to Semi-arid                 Cold to Cold 
##                          340                        95966 
## Dry subhumid to Dry subhumid        Dry subhumid to Humid 
##                         4022                          288 
##    Dry subhumid to Semi-arid        Humid to Dry subhumid 
##                         2242                         2484 
##               Humid to Humid           Humid to Semi-arid 
##                        57318                          198 
##            Hyperarid to Arid       Hyperarid to Hyperarid 
##                          240                         6402 
##                     NA to NA            Semi-arid to Arid 
##                         5130                         1209 
##    Semi-arid to Dry subhumid       Semi-arid to Semi-arid 
##                          277                        12514
cmip6s$change.cat.mean <- factor(cmip6s$change.cat.mean, 
                                  levels = c("Arid to Hyperarid","Semi-arid to Arid", "Dry subhumid to Semi-arid", "Humid to Semi-arid", "Humid to Dry subhumid", "Cold to Humid","Hyperarid to arid", "Arid to Semi-arid", "Dry subhumid to Humid", "Semi-arid to Dry subhumid", "Semi-arid to Humid"))

colors_changes <- c("Hyperarid to Arid" = "#E7D7B6", 
                     "Arid to Semi-arid" = colorvec[7],
                     "Dry subhumid to Humid" = colorvec[8],
                     "Semi-arid to Dry subhumid" = colorvec[9],
                     "Semi-arid to Humid" = colorvec[10],
                     "Humid to Cold" = colorvec[11],
                     "Cold to Humid" = "#E7D7B6",
                     "Humid to Dry subhumid" = colorvec[5],
                     "Humid to Semi-arid" = colorvec[4],
                     "Dry subhumid to Semi-arid" = colorvec[3],
                     "Semi-arid to Arid" = colorvec[2],
                     "Arid to Hyperarid" = colorvec[1])

8.2 Changes towards dryer categories cat mean

8.2.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + 
    geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmean == 1 & cat.dryer == "dryer"), 
                aes(x=lon, y = lat,  fill = change.cat.mean))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

8.2.1.1 Future


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmean == 1 & cat.dryer == "dryer"), 
                              aes(x=lon, y = lat,  fill = change.cat.mean))+
      borders(colour = "grey60")+
      scale_fill_manual(values = colors_changes)+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggarrange(plotlist = map_list, ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

8.3 Changes towards wetter categories cat mean

8.3.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + geom_tile(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & cat.dryer == "wetter"), 
                            aes(x=lon, y = lat,  fill = change.cat.mean))+
    borders(colour = "grey60", ylim = c(-40,40))+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+ylim(-50,50)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

8.3.2 Future


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmean == 1 & cat.dryer == "wetter"), 
                              aes(x=lon, y = lat,  fill = change.cat.mean))+
      borders(colour = "grey60", ylim = c(-40,40))+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

8.4 Big map of SSP585 2070-2100


g <- ggplot() + 
  geom_tile(data = subset(cmip6s, period == "2070_2100" & model == "SSP585" & change.catmean == 1), aes(x=lon, y = lat,  fill = change.cat.mean))+
      borders(colour = "grey60")+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      labs(fill = "Change of aridity category\nfor the period 2070-2100\nand SSP 5-8.5,\ncompared to 1970-2000", title = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 3))
   
print(g)

9 Tables

9.1 Changes of aridity category using cat mean

tab <- cmip6s %>% group_by(model, period, Continent) %>% summarise(Cold = table(cat.AI)[1], Hyperarid = table(cat.AI)[2], Arid = table(cat.AI)[3], Semiarid = table(cat.AI)[4], Drysubhumid = table(cat.AI)[5], Humid = table(cat.AI)[6], total = n()) 
write.table(tab, "tab.catAI.txt")

tab.future <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(period, model, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(period + model ~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

kable(tab.future[order(tab.future$model, decreasing = F),]) %>% kable_styling(bootstrap_options = "bordered") %>%
  column_spec(c(7,10), italic = T, include_thead = T) %>% row_spec(2, bold = T) 
model period Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
1 historical 1850_1880 3.4 6.0 7.1 3.3 19.8 28.8 48.9 2.6
2 historical 1970_2000 3.3 6.0 7.0 3.3 19.6 27.8 49.9 2.6
3 historical 1985_2015 3.3 6.1 6.9 3.3 19.6 27.3 50.6 2.6
4 SSP245 2030_2060 3.5 6.3 7.8 3.3 20.9 28.8 47.7 2.6
7 SSP245 2070_2100 3.6 6.4 7.8 3.5 21.3 29.4 46.6 2.6
5 SSP370 2030_2060 3.2 6.7 7.6 3.0 20.5 28.8 48.1 2.6
8 SSP370 2070_2100 4.0 6.3 8.3 3.7 22.3 28.5 46.6 2.6
6 SSP585 2030_2060 3.6 6.2 7.9 3.5 21.2 29.3 47.0 2.6
9 SSP585 2070_2100 3.8 6.4 8.4 3.6 22.2 30.1 45.1 2.6

# ggplot(tab.future, aes(x = period))+
#   geom_point(aes(y = Hyperarid, col = "Hyperarid", shape = model))+
#   geom_line(aes(y=Hyperarid, col = "Hyperarid", group = model, lty = model))+
#   geom_point(aes(y = Arid, col = "Arid", shape = model))+
#   geom_line(aes(y=Arid, col = "Arid", group = model, lty = model))+
#   geom_point(aes(y = get('Semi-arid'), col = "Semi-arid", shape = model))+
#   geom_line(aes(y=get('Semi-arid'), col = "Semi-arid", group = model, lty = model))+
#   geom_point(aes(y = get('Dry subhumid'), col = "Dry subhumid", shape = model))+
#   geom_line(aes(y=get('Dry subhumid'), col = "Dry subhumid", group = model, lty = model))+
#   geom_point(aes(y = Humid, col = "Humid", shape = model))+
#   geom_line(aes(y= Humid, col = "Humid", group = model, lty = model))+
#   scale_color_manual(values = col.cat)+
#   theme_minimal()

cmip6$cat.AI <- factor(cmip6$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))
tab.percent <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() %>% 
  group_by(period, model, cat.AI) %>% 
  summarise(mmmean = mean(count, na.rm = T), mmsd = sd(count, na.rm = T)) %>%
  mutate(percent = round(mmmean/sum(mmmean)*100, 1), sd.percent = round(mmsd/sum(mmmean)*100,1)) 
write.table(tab.percent, "tab.percent.txt")
tab.percent <- read.table("tab.percent.txt")

df245 <- rbind(subset(tab.percent, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.percent, model == "SSP245")) %>%
  group_by(period, model) %>% 
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

df245$cat.AI <- factor(df245$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

b245 <- ggplot(data = df245, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = paste(percent, " ± ", sd.percent, sep = "")), size = 3)+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 2-4.5", y = "%", x = "", col = "", fill = "")+
  theme_minimal()

df370 <- rbind(subset(tab.percent, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.percent, model == "SSP370")) %>%
  group_by(period, model) %>% 
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

df370$cat.AI <- factor(df370$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

b370 <- ggplot(data = df370, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = paste(percent, " ± ", sd.percent, sep = "")), size = 3)+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 3-7.0", y = "%", x = "", col = "", fill = "")+
  theme_minimal()


df585 <- rbind(subset(tab.percent, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.percent, model == "SSP585")) %>%
  group_by(period, model) %>% 
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

df585$cat.AI <- factor(df585$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

b585 <- ggplot(data = df585, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = paste(percent, " ± ", sd.percent, sep = "")), size = 3)+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 5-8.5", y = "%", x = "", col = "", fill = "")+
  theme_minimal()


ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)

9.2 Changes of aridity category using cat maj

tab.maj <- cmip6s %>% group_by(model, period, Continent) %>% summarise(Cold = table(cat.maj)[1], Hyperarid = table(cat.maj)[2], Arid = table(cat.maj)[3], Semiarid = table(cat.maj)[4], Drysubhumid = table(cat.maj)[5], Humid = table(cat.maj)[6], total = n()) 
write.table(tab.maj, "tab.catmaj.txt")

tab.future.maj <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(period, model, cat.maj) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(period + model ~cat.maj) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

kable(tab.future.maj[order(tab.future.maj$model, decreasing = F),]) %>% kable_styling(bootstrap_options = "bordered") %>%
  column_spec(c(7,10), italic = T, include_thead = T) %>% row_spec(2, bold = T) 
model period Hyperarid Arid Semi-arid Dry subhumid Sum drylands Humid Cold NA
1 historical 1850_1880 4.2 7.1 8.7 0.8 20.8 28.3 48.3 2.6
2 historical 1970_2000 4.1 7.1 8.6 0.7 20.5 28.1 48.9 2.6
3 historical 1985_2015 4.0 7.0 8.4 0.9 20.3 27.3 49.9 2.6
4 SSP245 2030_2060 4.0 7.8 9.3 0.8 21.9 28.4 47.2 2.6
7 SSP245 2070_2100 4.2 7.7 9.4 0.9 22.2 29.1 46.1 2.6
5 SSP370 2030_2060 4.1 7.2 9.4 0.8 21.5 28.6 47.3 2.6
8 SSP370 2070_2100 4.5 7.7 10.0 0.9 23.1 28.4 45.9 2.6
6 SSP585 2030_2060 4.2 7.5 9.6 1.0 22.3 28.8 46.3 2.6
9 SSP585 2070_2100 4.5 7.9 9.4 0.9 22.7 30.0 44.6 2.6
tab.flow.catmaj <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, cat.maj) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL)

df245 <- rbind(subset(tab.flow.catmaj, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow.catmaj, model == "SSP245")) %>%
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b245 <- ggplot(data = df245, aes(x = period, y = percent))+
  geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 2-4.5", y = "%", x = "")+
  theme_minimal()

df370 <- rbind(subset(tab.flow.catmaj, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow.catmaj, model == "SSP370")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b370 <- ggplot(data = df370, aes(x = period, y = percent))+
  geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 3-7.0", y = "%", x = "")+
  theme_minimal()


df585 <- rbind(subset(tab.flow.catmaj, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow.catmaj, model == "SSP585")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b585 <- ggplot(data = df585, aes(x = period, y = percent))+
  geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 5-8.5", y = "%", x = "")+
  theme_minimal()


ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)

9.3 Changes by continent, cat mean


cmip6$cat.AI <- factor(cmip6$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

tab.percent.cont <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, Continent, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() %>% 
  group_by(period, model, Continent, cat.AI) %>% 
  summarise(mmmean = mean(count, na.rm = T), mmsd = sd(count, na.rm = T)) %>%
  mutate(percent = round(mmmean/sum(mmmean)*100, 1), sd.percent = round(mmsd/sum(mmmean)*100,1)) 
write.table(tab.percent.cont, "tab.percent.continent.txt")
tab.percent.cont <- read.table("tab.percent.continent.txt")

tab_list_percent <- list()
tab_list_sd <- list()

for(i in unique(tab.percent.cont$Continent)){
  tab.i <- subset(tab.percent.cont, Continent == i)
  
  df.percent <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "percent")
  
  df.sd <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "sd.percent")
  # sd for a sum: square root of sum of squared sd
  
  drycats <- which(names(df.percent) %in% c("Hyperarid","Arid","Semi-arid","Dry subhumid"))
  
  df.percent <- df.percent %>% mutate("Sum drylands" = rowSums(.[drycats], na.rm = T))
  df.sd <- df.sd %>% mutate("Sum drylands" = round(sqrt(rowSums((.[drycats])^2, na.rm = T)),2))
  
  missing <- setdiff(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"), names(df.percent))
  
  df.percent[, missing] <- 0
  df.sd[,missing] <- 0
  
  df.percent <- df.percent %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.percent[is.na(df.percent)] <- 0
  
  tab_list_percent[[i]] <- df.percent
  
  df.sd <- df.sd %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.sd[is.na(df.sd)] <- 0
  
  tab_list_sd[[i]] <- df.sd
}

k_list <- list()


kdiff_list <- list()

for(i in names(tab_list_percent)){
  
  df <- tab_list_percent[[i]] %>% select(c("model", "period"))
  
  
  
  prop.ref <- tab_list_percent[[i]] %>% filter(period == "1970_2000")
  prop.ref <- prop.ref %>% rbind(prop.ref[rep(1,dim(df)[1]-1),])
  df2 <- tab_list_percent[[i]][,c(3:10)] - prop.ref[,c(3:10)] %>% round(1)
  diff.prop.percent <- cbind(df, df2)
  
  sd.ref <- tab_list_sd[[i]] %>% filter(period == "1970_2000")
  sd.ref <- sd.ref %>% rbind(sd.ref[rep(1,dim(df)[1]-1),])
  df3 <- sqrt(tab_list_sd[[i]][,c(3:10)]^2 + sd.ref[,c(3:10)]^2) %>% round(1)
  diff.sd.percent <- cbind(df, df3)
  
  df$Hyperarid <- paste(tab_list_percent[[i]]$Hyperarid, tab_list_sd[[i]]$Hyperarid, sep = " ± ")
  df$Arid <- paste(tab_list_percent[[i]]$Arid, tab_list_sd[[i]]$Arid, sep = " ± ")
  df$'Semi-Arid' <- paste(tab_list_percent[[i]][,5], tab_list_sd[[i]][,5], sep = " ± ")
  df$'Dry subhumid' <- paste(tab_list_percent[[i]][,6], tab_list_sd[[i]][,6], sep = " ± ")
  df$'Sum drylands' <- paste(tab_list_percent[[i]][,7], tab_list_sd[[i]][,7], sep = " ± ")
  df$Humid <- paste(tab_list_percent[[i]]$Humid, tab_list_sd[[i]]$Humid, sep = " ± ")
  df$Cold <- paste(tab_list_percent[[i]]$Cold, tab_list_sd[[i]]$Cold, sep = " ± ")
  
  df.diff <- df
  df.diff$Hyperarid <- paste(round(diff.prop.percent$Hyperarid, 1), round(diff.sd.percent$Hyperarid, 1), sep = " ± ")
  df.diff$Arid <- paste(round(diff.prop.percent$Arid, 1), round(diff.sd.percent$Arid, 1), sep = " ± ")
  df.diff$'Semi-Arid' <- paste(round(diff.prop.percent[,5], 1), round(diff.sd.percent[,5], 1), sep = " ± ")
  df.diff$'Dry subhumid' <- paste(round(diff.prop.percent[,6], 1), round(diff.sd.percent[,6], 1), sep = " ± ")
  df.diff$'Sum drylands' <- paste(round(diff.prop.percent[,7],1), round(diff.sd.percent[,7], 1), sep = " ± ")
  df.diff$Humid <- paste(round(diff.prop.percent$Humid, 1), round(diff.sd.percent$Humid, 1), sep = " ± ")
  df.diff$Cold <- paste(round(diff.prop.percent$Cold, 1), round(diff.sd.percent$Cold, 1), sep = " ± ")
  
  k_list[[i]] <- kable(df[order(df$model, decreasing = F),], caption = i) %>% kable_styling(bootstrap_options = "bordered") %>%
    column_spec(8, italic = T, background = colorvec[5]) %>% row_spec(2, bold = T) 
  kdiff_list[[i]] <- kable(df.diff[order(df.diff$model, decreasing = F),], caption = i) %>% kable_styling(bootstrap_options = "bordered") %>% 
    column_spec(8, italic = T, background = colorvec[5]) %>% row_spec(2, bold = T) 
  
}

9.3.1 Proportion of aridity categories by continent

k_list
$AFRICA
AFRICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 27.2 ± 4.7 17.1 ± 4.9 15.6 ± 3.6 6.3 ± 1.4 66.2 ± 7.81 33.7 ± 7.1 0 ± 0
2 historical 1970_2000 26.9 ± 4.7 16.7 ± 4.2 15.1 ± 3.4 6.2 ± 1.5 64.9 ± 7.32 32.8 ± 6.2 2.3 ± 0
3 historical 1985_2015 26.5 ± 4.8 17.1 ± 4.6 15.4 ± 3.3 6 ± 1.5 65 ± 7.57 32.7 ± 6.3 2.2 ± 0
4 SSP245 2030_2060 27 ± 4.8 18 ± 4.9 16.5 ± 3.5 6.4 ± 1.7 67.9 ± 7.89 32 ± 7.6 0 ± 0
7 SSP245 2070_2100 27.9 ± 4.7 17.4 ± 5 17.2 ± 3.4 6.2 ± 1.5 68.7 ± 7.8 31.3 ± 7.9 0 ± 0
5 SSP370 2030_2060 26 ± 5 16.6 ± 3.6 17.2 ± 1.8 6.8 ± 1 66.6 ± 6.5 32.3 ± 6.7 1.1 ± 0
8 SSP370 2070_2100 28 ± 3.5 17.6 ± 4.8 16.6 ± 3.1 6.2 ± 1.6 68.4 ± 6.89 30.4 ± 6.4 1.2 ± 0
6 SSP585 2030_2060 27.1 ± 4.6 17.6 ± 4.8 16.8 ± 3.5 6.3 ± 1.5 67.8 ± 7.66 32.1 ± 7.5 0 ± 0
9 SSP585 2070_2100 29.7 ± 10.5 17.3 ± 5.3 16.2 ± 4.2 6 ± 1.5 69.2 ± 12.58 30.9 ± 8.6 0 ± 0
$ASIA
ASIA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 3.9 ± 2.3 12.2 ± 1.9 9.2 ± 1.5 4 ± 1.4 29.3 ± 3.62 32.4 ± 7.7 38 ± 9.2
2 historical 1970_2000 3.8 ± 2.1 12 ± 1.8 8.9 ± 1.5 3.8 ± 1.1 28.5 ± 3.33 30.7 ± 8.5 40.5 ± 9.6
3 historical 1985_2015 3.7 ± 2.2 11.8 ± 1.6 8.6 ± 1.5 3.7 ± 1.1 27.8 ± 3.3 30.1 ± 7.6 41.8 ± 8.8
4 SSP245 2030_2060 3.5 ± 1.8 12.5 ± 1.6 9.4 ± 1.6 4.3 ± 1.7 29.7 ± 3.35 34.3 ± 8.1 35.8 ± 9.2
7 SSP245 2070_2100 3.7 ± 2 12.5 ± 1.7 9.5 ± 1.7 4.2 ± 1.7 29.9 ± 3.56 36.2 ± 8.6 33.7 ± 9.3
5 SSP370 2030_2060 3.4 ± 2.3 12.4 ± 1.5 8.9 ± 1.7 3.7 ± 1.3 28.4 ± 3.48 34 ± 7.2 37.3 ± 9.5
8 SSP370 2070_2100 4.6 ± 1.9 12.5 ± 1.5 9.6 ± 1.5 3.8 ± 1.2 30.5 ± 3.09 35.2 ± 9.2 34.1 ± 9.7
6 SSP585 2030_2060 3.8 ± 2 12.5 ± 1.6 9.5 ± 1.8 4.2 ± 1.6 30 ± 3.52 35.4 ± 8.4 34.4 ± 9.3
9 SSP585 2070_2100 4.8 ± 3.2 13.1 ± 3.2 10.4 ± 3.4 4.2 ± 1.5 32.5 ± 5.86 36.5 ± 11.7 30.7 ± 9.9
$CENTRAL-AMERICA
CENTRAL-AMERICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0 8.4 ± 3.6 30.1 ± 5.1 16.1 ± 8.7 54.6 ± 10.71 45 ± 6.1 0 ± 0
2 historical 1970_2000 0 ± 0 8.4 ± 4.1 29.7 ± 7.6 14.7 ± 7.8 52.8 ± 11.64 44.9 ± 12 2 ± 0
3 historical 1985_2015 0 ± 0 9 ± 4.2 30.7 ± 7 14.8 ± 7.4 54.5 ± 11.02 43.8 ± 11.5 1.4 ± 0
4 SSP245 2030_2060 0 ± 0 11.9 ± 6.2 36.3 ± 8 13.9 ± 6.2 62.1 ± 11.87 37.6 ± 8.6 0 ± 0
7 SSP245 2070_2100 0 ± 0 12.3 ± 5.5 38.1 ± 8.5 14.4 ± 7.9 64.8 ± 12.84 34.8 ± 8.6 0 ± 0
5 SSP370 2030_2060 0 ± 0 12.6 ± 6.3 31.4 ± 7.1 15 ± 7.9 59 ± 12.35 40.7 ± 11.8 0 ± 0
8 SSP370 2070_2100 4 ± 5.2 16 ± 5.5 35.5 ± 9.5 13.2 ± 5.2 68.7 ± 13.21 30.9 ± 10.6 0 ± 0
6 SSP585 2030_2060 0 ± 0 13 ± 7 37.6 ± 9 14.1 ± 7.1 64.7 ± 13.43 35 ± 9.4 0 ± 0
9 SSP585 2070_2100 2 ± 2.4 16.6 ± 8 41.2 ± 7.7 13.3 ± 7.7 73.1 ± 13.72 26.6 ± 8.8 0 ± 0
$EUROPE
EUROPE
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0 1.9 ± 1.2 9.6 ± 4.1 4.5 ± 2.8 16 ± 5.11 49 ± 15.3 34.8 ± 20.2
2 historical 1970_2000 0 ± 0 1.5 ± 1 8.9 ± 3.4 3.9 ± 2.1 14.3 ± 4.12 44.9 ± 16.1 40.6 ± 20.4
3 historical 1985_2015 0 ± 0 1.3 ± 1 8.4 ± 2.9 3.8 ± 2 13.5 ± 3.66 46 ± 16.1 40.3 ± 19.9
4 SSP245 2030_2060 0 ± 0 1.9 ± 0.9 10.3 ± 3.6 5.1 ± 1.9 17.3 ± 4.17 52.7 ± 15.4 29.8 ± 19.7
7 SSP245 2070_2100 0 ± 0 2.2 ± 1.1 10.1 ± 4.8 5.1 ± 2.3 17.4 ± 5.44 55.4 ± 14.4 27 ± 18.6
5 SSP370 2030_2060 0 ± 0 2.1 ± 1 9.7 ± 4.7 4.8 ± 2.5 16.6 ± 5.42 52.7 ± 15.7 30.4 ± 20
8 SSP370 2070_2100 0 ± 0 3.1 ± 1.1 11.4 ± 5.9 6.4 ± 3.8 20.9 ± 7.1 52.8 ± 14.3 26 ± 18.6
6 SSP585 2030_2060 0 ± 0 2.3 ± 1.1 10.5 ± 3.7 4.9 ± 2.4 17.7 ± 4.55 55.2 ± 15.8 26.9 ± 19.9
9 SSP585 2070_2100 0 ± 0 3 ± 1.7 11.1 ± 5.4 6.1 ± 2.9 20.2 ± 6.36 57.2 ± 12.8 22.5 ± 17.3
$EUROPE-AFRICA
EUROPE-AFRICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 5.9 ± 4.9 23.6 ± 7 21.3 ± 8.8 8.8 ± 3.3 59.6 ± 12.7 24.9 ± 14.2 14.5 ± 0
2 historical 1970_2000 7 ± 5.5 26.5 ± 7.5 23.3 ± 10.3 9.9 ± 3.2 66.7 ± 14.24 27.5 ± 16.5 4.5 ± 7.5
3 historical 1985_2015 7.1 ± 5.9 26.3 ± 7.1 22.2 ± 7.7 10.3 ± 4.2 65.9 ± 12.73 27.8 ± 14.3 5.2 ± 7
4 SSP245 2030_2060 9.9 ± 6.6 28.9 ± 8.6 26.7 ± 8.9 9.4 ± 3.4 74.9 ± 14.43 23.1 ± 16.2 0.9 ± 0
7 SSP245 2070_2100 11.7 ± 5.7 28.6 ± 8.8 27.2 ± 8.1 9.5 ± 3.1 77 ± 13.61 21.9 ± 13.9 0 ± 0
5 SSP370 2030_2060 10.3 ± 6.3 28.2 ± 8.4 26 ± 8.7 9 ± 2.8 73.5 ± 13.92 24.8 ± 14 0.5 ± 0
8 SSP370 2070_2100 13.2 ± 6.2 27.3 ± 8.6 28.3 ± 9.7 8.7 ± 3.6 77.5 ± 14.81 21.2 ± 16.2 0 ± 0
6 SSP585 2030_2060 11.1 ± 6.4 28.8 ± 8.3 27.5 ± 8.8 8.7 ± 3 76.1 ± 14.01 22.6 ± 14.1 0.2 ± 0
9 SSP585 2070_2100 15.2 ± 5.7 29.1 ± 10.1 27.9 ± 8.8 8 ± 2.8 80.2 ± 14.82 18.6 ± 13.6 0 ± 0
$NORTH-AMERICA
NORTH-AMERICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0 0.4 ± 0.2 5.8 ± 5.6 4.3 ± 2.6 10.5 ± 6.18 35.1 ± 9 54.2 ± 10.9
2 historical 1970_2000 0 ± 0 0.4 ± 0.2 4.3 ± 4.4 4.2 ± 2.4 8.9 ± 5.02 34.5 ± 9.8 56.5 ± 11.7
3 historical 1985_2015 0 ± 0 0.3 ± 0.2 4.6 ± 3.8 4.1 ± 2.2 9 ± 4.4 33.6 ± 10.3 57.2 ± 11.5
4 SSP245 2030_2060 0 ± 0 0.4 ± 0.3 6.5 ± 5.5 4.7 ± 2.6 11.6 ± 6.09 37.5 ± 8.7 50.8 ± 11.1
7 SSP245 2070_2100 0 ± 0 0.5 ± 0.3 6.4 ± 5.1 5.1 ± 2.9 12 ± 5.87 39.7 ± 9 48.1 ± 11
5 SSP370 2030_2060 0 ± 0 0.5 ± 0.4 5.9 ± 4.8 4.6 ± 2.7 11 ± 5.52 37 ± 9.7 51.8 ± 12.2
8 SSP370 2070_2100 0.1 ± 0 0.8 ± 0.7 7.6 ± 4 5.2 ± 2 13.7 ± 4.53 38.4 ± 10.6 47.7 ± 11.2
6 SSP585 2030_2060 0 ± 0 0.5 ± 0.4 6.6 ± 5.2 5.2 ± 2.8 12.3 ± 5.92 38.2 ± 9.1 49.3 ± 11.7
9 SSP585 2070_2100 0.2 ± 0 1 ± 1.4 7 ± 6.1 5.2 ± 2.5 13.4 ± 6.74 42.1 ± 9.3 44.5 ± 11.3
$OCEANIA
OCEANIA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 5.3 ± 0 32 ± 26 36.7 ± 17.1 6.6 ± 3.1 80.6 ± 31.27 18.9 ± 12.5 0.3 ± 0.2
2 historical 1970_2000 5.2 ± 0 31.3 ± 26.4 36.8 ± 17.4 6.8 ± 3.1 80.1 ± 31.77 19.3 ± 12.8 0.5 ± 0.4
3 historical 1985_2015 5.2 ± 0 33.4 ± 26.1 35.5 ± 17.9 6.8 ± 3.2 80.9 ± 31.81 18.4 ± 12.3 0.6 ± 0.3
4 SSP245 2030_2060 10.4 ± 0 35.8 ± 23.5 31.7 ± 15.1 5.7 ± 2.6 83.6 ± 28.05 15.9 ± 11.9 0.3 ± 0.3
7 SSP245 2070_2100 10.5 ± 0 37.9 ± 22.2 29.1 ± 13.1 6.2 ± 3.3 83.7 ± 25.99 16.2 ± 12 0.1 ± 0
5 SSP370 2030_2060 10.7 ± 0 34.8 ± 22.6 33.4 ± 15.3 5.8 ± 2.5 84.7 ± 27.41 14.9 ± 8.1 0.3 ± 0.2
8 SSP370 2070_2100 10.5 ± 10.2 38.1 ± 19.7 28.5 ± 13.8 5.7 ± 3.2 82.8 ± 26.32 16.8 ± 16.3 0.3 ± 0.3
6 SSP585 2030_2060 5.6 ± 3.6 35.6 ± 25.3 33.8 ± 16.7 6.7 ± 3.6 81.7 ± 30.74 17.7 ± 13.8 0.3 ± 0.3
9 SSP585 2070_2100 9.3 ± 6.4 39.4 ± 21.7 29.4 ± 12.9 5.5 ± 2.5 83.6 ± 26.16 16.2 ± 12.5 0 ± 0
$SOUTH-AMERICA
SOUTH-AMERICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0.1 ± 0.1 2.1 ± 1.7 15.2 ± 8.7 7.4 ± 2.6 24.8 ± 9.24 73.2 ± 11.7 1.9 ± 1
2 historical 1970_2000 0.1 ± 0.1 1.8 ± 1.3 13.9 ± 4.3 7.3 ± 2.3 23.1 ± 5.05 74.5 ± 7.3 2.3 ± 2.6
3 historical 1985_2015 0.2 ± 0.1 2 ± 1.2 14.5 ± 5.1 7.5 ± 2.4 24.2 ± 5.76 73.6 ± 7.7 2.1 ± 2.6
4 SSP245 2030_2060 0.2 ± 0.2 2.7 ± 2.2 17.4 ± 9.6 8.8 ± 3.6 29.1 ± 10.49 69.4 ± 13.1 1.3 ± 0.9
7 SSP245 2070_2100 0.2 ± 0.1 3.3 ± 3.5 18.3 ± 9.1 9.3 ± 4.6 31.1 ± 10.78 67.5 ± 14.7 1.3 ± 0.7
5 SSP370 2030_2060 0.2 ± 0.1 2.2 ± 1.7 16.3 ± 5.5 8.8 ± 3.7 27.5 ± 6.84 70.5 ± 8.9 1.9 ± 1.6
8 SSP370 2070_2100 0.2 ± 0.1 3.3 ± 2.5 18.6 ± 7 8.7 ± 3.8 30.8 ± 8.35 67.3 ± 10.8 1.8 ± 1.8
6 SSP585 2030_2060 0.2 ± 0.1 3.1 ± 2.5 18.6 ± 9 9.2 ± 4.4 31.1 ± 10.33 67.6 ± 13.1 1.2 ± 0.8
9 SSP585 2070_2100 0.6 ± 1 4.4 ± 3.8 19.6 ± 9.5 9.7 ± 3.9 34.3 ± 11 64.4 ± 12.6 1.1 ± 0.7

9.3.2 Changes in proportion of aridity categories by continent

kdiff_list
$AFRICA
AFRICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0.3 ± 6.6 0.4 ± 6.5 0.5 ± 5 0.1 ± 2.1 1.3 ± 10.7 0.9 ± 9.4 -2.3 ± 0
2 historical 1970_2000 0 ± 6.6 0 ± 5.9 0 ± 4.8 0 ± 2.1 0 ± 10.4 0 ± 8.8 0 ± 0
3 historical 1985_2015 -0.4 ± 6.7 0.4 ± 6.2 0.3 ± 4.7 -0.2 ± 2.1 0.1 ± 10.5 -0.1 ± 8.8 -0.1 ± 0
4 SSP245 2030_2060 0.1 ± 6.7 1.3 ± 6.5 1.4 ± 4.9 0.2 ± 2.3 3 ± 10.8 -0.8 ± 9.8 -2.3 ± 0
7 SSP245 2070_2100 1 ± 6.6 0.7 ± 6.5 2.1 ± 4.8 0 ± 2.1 3.8 ± 10.7 -1.5 ± 10 -2.3 ± 0
5 SSP370 2030_2060 -0.9 ± 6.9 -0.1 ± 5.5 2.1 ± 3.8 0.6 ± 1.8 1.7 ± 9.8 -0.5 ± 9.1 -1.2 ± 0
8 SSP370 2070_2100 1.1 ± 5.9 0.9 ± 6.4 1.5 ± 4.6 0 ± 2.2 3.5 ± 10.1 -2.4 ± 8.9 -1.1 ± 0
6 SSP585 2030_2060 0.2 ± 6.6 0.9 ± 6.4 1.7 ± 4.9 0.1 ± 2.1 2.9 ± 10.6 -0.7 ± 9.7 -2.3 ± 0
9 SSP585 2070_2100 2.8 ± 11.5 0.6 ± 6.8 1.1 ± 5.4 -0.2 ± 2.1 4.3 ± 14.6 -1.9 ± 10.6 -2.3 ± 0
$ASIA
ASIA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0.1 ± 3.1 0.2 ± 2.6 0.3 ± 2.1 0.2 ± 1.8 0.8 ± 4.9 1.7 ± 11.5 -2.5 ± 13.3
2 historical 1970_2000 0 ± 3 0 ± 2.5 0 ± 2.1 0 ± 1.6 0 ± 4.7 0 ± 12 0 ± 13.6
3 historical 1985_2015 -0.1 ± 3 -0.2 ± 2.4 -0.3 ± 2.1 -0.1 ± 1.6 -0.7 ± 4.7 -0.6 ± 11.4 1.3 ± 13
4 SSP245 2030_2060 -0.3 ± 2.8 0.5 ± 2.4 0.5 ± 2.2 0.5 ± 2 1.2 ± 4.7 3.6 ± 11.7 -4.7 ± 13.3
7 SSP245 2070_2100 -0.1 ± 2.9 0.5 ± 2.5 0.6 ± 2.3 0.4 ± 2 1.4 ± 4.9 5.5 ± 12.1 -6.8 ± 13.4
5 SSP370 2030_2060 -0.4 ± 3.1 0.4 ± 2.3 0 ± 2.3 -0.1 ± 1.7 -0.1 ± 4.8 3.3 ± 11.1 -3.2 ± 13.5
8 SSP370 2070_2100 0.8 ± 2.8 0.5 ± 2.3 0.7 ± 2.1 0 ± 1.6 2 ± 4.5 4.5 ± 12.5 -6.4 ± 13.6
6 SSP585 2030_2060 0 ± 2.9 0.5 ± 2.4 0.6 ± 2.3 0.4 ± 1.9 1.5 ± 4.8 4.7 ± 12 -6.1 ± 13.4
9 SSP585 2070_2100 1 ± 3.8 1.1 ± 3.7 1.5 ± 3.7 0.4 ± 1.9 4 ± 6.7 5.8 ± 14.5 -9.8 ± 13.8
$CENTRAL-AMERICA
CENTRAL-AMERICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0 0 ± 5.5 0.4 ± 9.2 1.4 ± 11.7 1.8 ± 15.8 0.1 ± 13.5 -2 ± 0
2 historical 1970_2000 0 ± 0 0 ± 5.8 0 ± 10.7 0 ± 11 0 ± 16.5 0 ± 17 0 ± 0
3 historical 1985_2015 0 ± 0 0.6 ± 5.9 1 ± 10.3 0.1 ± 10.8 1.7 ± 16 -1.1 ± 16.6 -0.6 ± 0
4 SSP245 2030_2060 0 ± 0 3.5 ± 7.4 6.6 ± 11 -0.8 ± 10 9.3 ± 16.6 -7.3 ± 14.8 -2 ± 0
7 SSP245 2070_2100 0 ± 0 3.9 ± 6.9 8.4 ± 11.4 -0.3 ± 11.1 12 ± 17.3 -10.1 ± 14.8 -2 ± 0
5 SSP370 2030_2060 0 ± 0 4.2 ± 7.5 1.7 ± 10.4 0.3 ± 11.1 6.2 ± 17 -4.2 ± 16.8 -2 ± 0
8 SSP370 2070_2100 4 ± 5.2 7.6 ± 6.9 5.8 ± 12.2 -1.5 ± 9.4 15.9 ± 17.6 -14 ± 16 -2 ± 0
6 SSP585 2030_2060 0 ± 0 4.6 ± 8.1 7.9 ± 11.8 -0.6 ± 10.5 11.9 ± 17.8 -9.9 ± 15.2 -2 ± 0
9 SSP585 2070_2100 2 ± 2.4 8.2 ± 9 11.5 ± 10.8 -1.4 ± 11 20.3 ± 18 -18.3 ± 14.9 -2 ± 0
$EUROPE
EUROPE
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0 0.4 ± 1.6 0.7 ± 5.3 0.6 ± 3.5 1.7 ± 6.6 4.1 ± 22.2 -5.8 ± 28.7
2 historical 1970_2000 0 ± 0 0 ± 1.4 0 ± 4.8 0 ± 3 0 ± 5.8 0 ± 22.8 0 ± 28.8
3 historical 1985_2015 0 ± 0 -0.2 ± 1.4 -0.5 ± 4.5 -0.1 ± 2.9 -0.8 ± 5.5 1.1 ± 22.8 -0.3 ± 28.5
4 SSP245 2030_2060 0 ± 0 0.4 ± 1.3 1.4 ± 5 1.2 ± 2.8 3 ± 5.9 7.8 ± 22.3 -10.8 ± 28.4
7 SSP245 2070_2100 0 ± 0 0.7 ± 1.5 1.2 ± 5.9 1.2 ± 3.1 3.1 ± 6.8 10.5 ± 21.6 -13.6 ± 27.6
5 SSP370 2030_2060 0 ± 0 0.6 ± 1.4 0.8 ± 5.8 0.9 ± 3.3 2.3 ± 6.8 7.8 ± 22.5 -10.2 ± 28.6
8 SSP370 2070_2100 0 ± 0 1.6 ± 1.5 2.5 ± 6.8 2.5 ± 4.3 6.6 ± 8.2 7.9 ± 21.5 -14.6 ± 27.6
6 SSP585 2030_2060 0 ± 0 0.8 ± 1.5 1.6 ± 5 1 ± 3.2 3.4 ± 6.1 10.3 ± 22.6 -13.7 ± 28.5
9 SSP585 2070_2100 0 ± 0 1.5 ± 2 2.2 ± 6.4 2.2 ± 3.6 5.9 ± 7.6 12.3 ± 20.6 -18.1 ± 26.7
$EUROPE-AFRICA
EUROPE-AFRICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 -1.1 ± 7.4 -2.9 ± 10.3 -2 ± 13.5 -1.1 ± 4.6 -7.1 ± 19.1 -2.6 ± 21.8 10 ± 7.5
2 historical 1970_2000 0 ± 7.8 0 ± 10.6 0 ± 14.6 0 ± 4.5 0 ± 20.1 0 ± 23.3 0 ± 10.6
3 historical 1985_2015 0.1 ± 8.1 -0.2 ± 10.3 -1.1 ± 12.9 0.4 ± 5.3 -0.8 ± 19.1 0.3 ± 21.8 0.7 ± 10.3
4 SSP245 2030_2060 2.9 ± 8.6 2.4 ± 11.4 3.4 ± 13.6 -0.5 ± 4.7 8.2 ± 20.3 -4.4 ± 23.1 -3.6 ± 7.5
7 SSP245 2070_2100 4.7 ± 7.9 2.1 ± 11.6 3.9 ± 13.1 -0.4 ± 4.5 10.3 ± 19.7 -5.6 ± 21.6 -4.5 ± 7.5
5 SSP370 2030_2060 3.3 ± 8.4 1.7 ± 11.3 2.7 ± 13.5 -0.9 ± 4.3 6.8 ± 19.9 -2.7 ± 21.6 -4 ± 7.5
8 SSP370 2070_2100 6.2 ± 8.3 0.8 ± 11.4 5 ± 14.1 -1.2 ± 4.8 10.8 ± 20.5 -6.3 ± 23.1 -4.5 ± 7.5
6 SSP585 2030_2060 4.1 ± 8.4 2.3 ± 11.2 4.2 ± 13.5 -1.2 ± 4.4 9.4 ± 20 -4.9 ± 21.7 -4.3 ± 7.5
9 SSP585 2070_2100 8.2 ± 7.9 2.6 ± 12.6 4.6 ± 13.5 -1.9 ± 4.3 13.5 ± 20.6 -8.9 ± 21.4 -4.5 ± 7.5
$NORTH-AMERICA
NORTH-AMERICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0 0 ± 0.3 1.5 ± 7.1 0.1 ± 3.5 1.6 ± 8 0.6 ± 13.3 -2.3 ± 16
2 historical 1970_2000 0 ± 0 0 ± 0.3 0 ± 6.2 0 ± 3.4 0 ± 7.1 0 ± 13.9 0 ± 16.5
3 historical 1985_2015 0 ± 0 -0.1 ± 0.3 0.3 ± 5.8 -0.1 ± 3.3 0.1 ± 6.7 -0.9 ± 14.2 0.7 ± 16.4
4 SSP245 2030_2060 0 ± 0 0 ± 0.4 2.2 ± 7 0.5 ± 3.5 2.7 ± 7.9 3 ± 13.1 -5.7 ± 16.1
7 SSP245 2070_2100 0 ± 0 0.1 ± 0.4 2.1 ± 6.7 0.9 ± 3.8 3.1 ± 7.7 5.2 ± 13.3 -8.4 ± 16.1
5 SSP370 2030_2060 0 ± 0 0.1 ± 0.4 1.6 ± 6.5 0.4 ± 3.6 2.1 ± 7.5 2.5 ± 13.8 -4.7 ± 16.9
8 SSP370 2070_2100 0.1 ± 0 0.4 ± 0.7 3.3 ± 5.9 1 ± 3.1 4.8 ± 6.8 3.9 ± 14.4 -8.8 ± 16.2
6 SSP585 2030_2060 0 ± 0 0.1 ± 0.4 2.3 ± 6.8 1 ± 3.7 3.4 ± 7.8 3.7 ± 13.4 -7.2 ± 16.5
9 SSP585 2070_2100 0.2 ± 0 0.6 ± 1.4 2.7 ± 7.5 1 ± 3.5 4.5 ± 8.4 7.6 ± 13.5 -12 ± 16.3
$OCEANIA
OCEANIA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0.1 ± 0 0.7 ± 37.1 -0.1 ± 24.4 -0.2 ± 4.4 0.5 ± 44.6 -0.4 ± 17.9 -0.2 ± 0.4
2 historical 1970_2000 0 ± 0 0 ± 37.3 0 ± 24.6 0 ± 4.4 0 ± 44.9 0 ± 18.1 0 ± 0.6
3 historical 1985_2015 0 ± 0 2.1 ± 37.1 -1.3 ± 25 0 ± 4.5 0.8 ± 45 -0.9 ± 17.8 0.1 ± 0.5
4 SSP245 2030_2060 5.2 ± 0 4.5 ± 35.3 -5.1 ± 23 -1.1 ± 4 3.5 ± 42.4 -3.4 ± 17.5 -0.2 ± 0.5
7 SSP245 2070_2100 5.3 ± 0 6.6 ± 34.5 -7.7 ± 21.8 -0.6 ± 4.5 3.6 ± 41 -3.1 ± 17.5 -0.4 ± 0.4
5 SSP370 2030_2060 5.5 ± 0 3.5 ± 34.8 -3.4 ± 23.2 -1 ± 4 4.6 ± 42 -4.4 ± 15.1 -0.2 ± 0.4
8 SSP370 2070_2100 5.3 ± 10.2 6.8 ± 32.9 -8.3 ± 22.2 -1.1 ± 4.5 2.7 ± 41.3 -2.5 ± 20.7 -0.2 ± 0.5
6 SSP585 2030_2060 0.4 ± 3.6 4.3 ± 36.6 -3 ± 24.1 -0.1 ± 4.8 1.6 ± 44.2 -1.6 ± 18.8 -0.2 ± 0.5
9 SSP585 2070_2100 4.1 ± 6.4 8.1 ± 34.2 -7.4 ± 21.7 -1.3 ± 4 3.5 ± 41.2 -3.1 ± 17.9 -0.5 ± 0.4
$SOUTH-AMERICA
SOUTH-AMERICA
model period Hyperarid Arid Semi-Arid Dry subhumid Sum drylands Humid Cold
1 historical 1850_1880 0 ± 0.1 0.3 ± 2.1 1.3 ± 9.7 0.1 ± 3.5 1.7 ± 10.5 -1.3 ± 13.8 -0.4 ± 2.8
2 historical 1970_2000 0 ± 0.1 0 ± 1.8 0 ± 6.1 0 ± 3.3 0 ± 7.1 0 ± 10.3 0 ± 3.7
3 historical 1985_2015 0.1 ± 0.1 0.2 ± 1.8 0.6 ± 6.7 0.2 ± 3.3 1.1 ± 7.7 -0.9 ± 10.6 -0.2 ± 3.7
4 SSP245 2030_2060 0.1 ± 0.2 0.9 ± 2.6 3.5 ± 10.5 1.5 ± 4.3 6 ± 11.6 -5.1 ± 15 -1 ± 2.8
7 SSP245 2070_2100 0.1 ± 0.1 1.5 ± 3.7 4.4 ± 10.1 2 ± 5.1 8 ± 11.9 -7 ± 16.4 -1 ± 2.7
5 SSP370 2030_2060 0.1 ± 0.1 0.4 ± 2.1 2.4 ± 7 1.5 ± 4.4 4.4 ± 8.5 -4 ± 11.5 -0.4 ± 3.1
8 SSP370 2070_2100 0.1 ± 0.1 1.5 ± 2.8 4.7 ± 8.2 1.4 ± 4.4 7.7 ± 9.8 -7.2 ± 13 -0.5 ± 3.2
6 SSP585 2030_2060 0.1 ± 0.1 1.3 ± 2.8 4.7 ± 10 1.9 ± 5 8 ± 11.5 -6.9 ± 15 -1.1 ± 2.7
9 SSP585 2070_2100 0.5 ± 1 2.6 ± 4 5.7 ± 10.4 2.4 ± 4.5 11.2 ± 12.1 -10.1 ± 14.6 -1.2 ± 2.7

9.3.3 Changes in proportion of aridity categories, by continent, calculated by model


cmip6$cat.AI <- factor(cmip6$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

tab.percent.cont <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, Continent, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() %>% 
  group_by(period, model, Continent, cat.AI) %>% 
  summarise(mmmean = mean(count, na.rm = T), mmsd = sd(count, na.rm = T)) %>%
  mutate(percent = round(mmmean/sum(mmmean)*100, 1), sd.percent = round(mmsd/sum(mmmean)*100,1)) 
write.table(tab.percent.cont, "tab.percent.continent.txt")

tab.cont <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, Continent, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() 

list_ref <- list()

for(i in unique(tab.cont$source)){
  tab.cont.i <- filter(tab.cont, source == i)
  tab.cont.i.ref <- filter(tab.cont.i, period == "1970_2000")
  
  for(j in unique(tab.cont.i$Continent)){
    tab.cont.i.j <- filter(tab.cont.i, Continent == j)
    tab.cont.i.j.ref <- filter(tab.cont.i.ref, Continent == j)
    
    for(k in unique(tab.cont.i$cat.AI)){
      
      index <- paste(i, j, k, sep = "_")
      tab.cont.i.j.k <- filter(tab.cont.i.j, cat.AI == k)
      
      lines <- ifelse(dim(tab.cont.i.j.k)[1] > 0, dim(tab.cont.i.j.k)[1] - 1, 1)
      
      tab.cont.i.j.k.ref <- subset(tab.cont.i.j.ref, cat.AI == k)
      
      if(dim(tab.cont.i.j.k.ref)[1] == 0){
        names1 <- names(tab.cont.i.j.k.ref)
        line1 <- c("1970_2000","historical", i, j, k, 0)
        tab.cont.i.j.k.ref <- rbind(tab.cont.i.j.k.ref, line1) %>%
          setNames(names1)
        
      }else{
        tab.cont.i.j.k.ref <- tab.cont.i.j.k.ref
      }
      
      tab.cont.i.j.k.ref <- tab.cont.i.j.k.ref %>% rbind(tab.cont.i.j.k.ref[rep(1,lines),])
      tab.cont.i.j.k$diff <- tab.cont.i.j.k$count - as.numeric(tab.cont.i.j.k.ref$count)  
      
      list_ref[[index]] <- tab.cont.i.j.k
    }
  }
}


tab.cont.diff <- bind_rows(list_ref, .id = "column_label") %>%
  group_by(period, model, Continent, cat.AI) %>% 
  summarise(diff.mean = mean(diff, na.rm = T), diff.sd = sd(diff, na.rm = T))
write.table(tab.cont.diff, "tab.cont.diff.txt")
tab.cont.diff <- read.table("tab.cont.diff.txt")

tab_list_percent <- list()
tab_list_sd <- list()

for(i in unique(tab.cont.diff$Continent)){
  tab.i <- subset(tab.cont.diff, Continent == i)
  
  df.percent <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "diff.mean")
  
  df.sd <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "diff.sd")
  # sd for a sum: square root of sum of squared sd
  
  drycats <- which(names(df.percent) %in% c("Hyperarid","Arid","Semi-arid","Dry subhumid"))
  
  df.percent <- df.percent %>% mutate("Sum drylands" = rowSums(.[drycats], na.rm = T))
  df.sd <- df.sd %>% mutate("Sum drylands" = sqrt(rowSums((.[drycats])^2, na.rm = T)),2)
  
  missing <- setdiff(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"), names(df.percent))
  
  df.percent[, missing] <- 0
  df.sd[,missing] <- 0
  
  df.percent <- df.percent %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.percent[is.na(df.percent)] <- 0
  df.percent[,c(3:9)] <- round(df.percent[,c(3:9)], 1)
  
  tab_list_percent[[i]] <- df.percent
  
  df.sd <- df.sd %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.sd[is.na(df.sd)] <- 0
  df.sd[ , c(3:9)] <- round(df.sd[, c(3:9)], 1)
  
  tab_list_sd[[i]] <- df.sd
}

k_list <- list()

for(i in names(tab_list_percent)){
  
  df <- tab_list_percent[[i]] %>% select(c("model", "period"))
  
  df$Hyperarid <- paste(tab_list_percent[[i]]$Hyperarid, tab_list_sd[[i]]$Hyperarid, sep = " ± ")
  df$Arid <- paste(tab_list_percent[[i]]$Arid, tab_list_sd[[i]]$Arid, sep = " ± ")
  df$'Semi-Arid' <- paste(tab_list_percent[[i]][,5], tab_list_sd[[i]][,5], sep = " ± ")
  df$'Dry subhumid' <- paste(tab_list_percent[[i]][,6], tab_list_sd[[i]][,6], sep = " ± ")
  df$'Sum drylands' <- paste(tab_list_percent[[i]][,7], tab_list_sd[[i]][,7], sep = " ± ")
  df$Humid <- paste(tab_list_percent[[i]]$Humid, tab_list_sd[[i]]$Humid, sep = " ± ")
  df$Cold <- paste(tab_list_percent[[i]]$Cold, tab_list_sd[[i]]$Cold, sep = " ± ")
  
  k_list[[i]] <- kable(df[order(df$model, decreasing = F)], caption = i, digits = 2) %>% kable_styling(bootstrap_options = "bordered") %>%
    column_spec(8, italic = T, background = colorvec[5]) %>% row_spec(2, bold = T) 
  
}

k_list
$AFRICA
AFRICA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 -5.8 ± 24.1 3 ± 22.7 1.5 ± 48.9 4.8 ± 34 2.8 ± 38.7 -0.4 ± 12 0 ± 0
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 -10.5 ± 19.5 11 ± 21.6 2.8 ± 33.1 5.9 ± 11.9 -2.7 ± 15.5 -3.5 ± 10.1 -2 ± 0
SSP245 2030_2060 -10.5 ± 41 24.3 ± 29.7 42.2 ± 66.6 25.8 ± 40.3 -37.8 ± 53.2 2.5 ± 15.5 0 ± 0
SSP370 2030_2060 -29.5 ± 88.4 -5.2 ± 60.9 23.7 ± 145 45.5 ± 86.1 -21.4 ± 43.4 12.9 ± 45.6 -30 ± 0
SSP585 2030_2060 -7.8 ± 50.9 13.6 ± 39.6 39 ± 82.2 33.2 ± 47.6 -34.6 ± 52.7 0 ± 18.4 0 ± 0
SSP245 2070_2100 10.1 ± 66.8 9 ± 46.2 58.8 ± 99 42.1 ± 51.9 -54.5 ± 61.7 -2.3 ± 22.7 0 ± 0
SSP370 2070_2100 19.5 ± 66.3 18.5 ± 58.2 69.2 ± 103.5 31.1 ± 46.6 -67.1 ± 52.3 0.2 ± 27.4 -28 ± 0
SSP585 2070_2100 45.6 ± 209.6 9.2 ± 79.8 67.5 ± 240.3 16.5 ± 83.8 -62.8 ± 98 -3.8 ± 20.8 0 ± 0
$ASIA
ASIA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 4.3 ± 28 8.8 ± 41.5 40.2 ± 88.6 11.3 ± 61.5 90.2 ± 176.8 15.8 ± 39.6 -130.5 ± 287.7
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 -5.2 ± 20 -12.6 ± 42 -36.4 ± 57.8 -17 ± 29.8 -33.6 ± 169.1 -1.6 ± 17.2 70 ± 234.1
SSP245 2030_2060 -15.7 ± 38.1 24.2 ± 40.5 57.8 ± 123.7 22.2 ± 88.9 192.7 ± 202.9 27 ± 65.5 -250.5 ± 375.8
SSP370 2030_2060 -23.8 ± 61.4 21.6 ± 45.7 -4.5 ± 98.3 0 ± 55.9 174.2 ± 120.1 -2.3 ± 26 -169.7 ± 100.4
SSP585 2030_2060 -1.8 ± 30.8 23.8 ± 34.1 72.9 ± 115.9 27.1 ± 91.2 250 ± 201.3 23.8 ± 54.9 -322.9 ± 358.7
SSP245 2070_2100 -8.5 ± 31.9 23.6 ± 37.9 68.1 ± 133.8 30.1 ± 101.2 294.5 ± 208.7 22.8 ± 72.2 -362.5 ± 381.8
SSP370 2070_2100 40.1 ± 64.3 25.3 ± 63.8 101.2 ± 108.9 32.7 ± 53 241.2 ± 179.9 3.1 ± 29 -342.3 ± 156
SSP585 2070_2100 55.1 ± 178.4 59.5 ± 146 215 ± 313.2 77.9 ± 202.8 324.1 ± 262.6 22.5 ± 61.9 -539.1 ± 416.4
$CENTRAL-AMERICA
CENTRAL-AMERICA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 0 ± 0 -0.3 ± 4.7 2.5 ± 21.5 -0.5 ± 18 -2 ± 28.2 3.3 ± 10.8 0 ± 0
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 0 ± 0 1.7 ± 4 4.2 ± 8.6 2.5 ± 6.3 -4 ± 5 -0.1 ± 4.3 -2 ± 0
SSP245 2030_2060 0 ± 0 9.8 ± 11.1 24 ± 29.7 17.3 ± 25.7 -23.5 ± 26.8 -3.2 ± 10 0 ± 0
SSP370 2030_2060 0 ± 0 11.8 ± 12.4 15.1 ± 31.2 3.3 ± 24 -14.6 ± 37.3 0 ± 15.6 0 ± 0
SSP585 2030_2060 0 ± 0 12.8 ± 14 31.5 ± 36.2 21.1 ± 29.1 -31 ± 25.5 -2.5 ± 16.4 0 ± 0
SSP245 2070_2100 0 ± 0 11 ± 10.3 31.9 ± 36.4 22.6 ± 29.2 -31.5 ± 27.3 -1.7 ± 19 0 ± 0
SSP370 2070_2100 12 ± 15.6 23.2 ± 15.1 50.2 ± 39 18.7 ± 29.5 -39.6 ± 23.8 -3.7 ± 13.4 0 ± 0
SSP585 2070_2100 6 ± 7.1 25.2 ± 24.9 60.6 ± 42.8 35.6 ± 28.6 -55.1 ± 42.4 -6.2 ± 18.6 0 ± 0
$EUROPE
EUROPE
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 0 ± 0 4.7 ± 11.2 20.7 ± 20.7 8.2 ± 13.3 50.3 ± 106.6 7.8 ± 11.3 -70 ± 115
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 0 ± 0 -2.2 ± 9.8 -9.2 ± 27.8 -5.8 ± 24.6 12.9 ± 48.5 -1.2 ± 8.5 -4.3 ± 73.1
SSP245 2030_2060 0 ± 0 4.7 ± 7.4 33.7 ± 18.9 17.5 ± 14.6 97.5 ± 126 11.5 ± 9.6 -128.5 ± 146.3
SSP370 2030_2060 0 ± 0 7.2 ± 8.5 35.2 ± 16.1 17.6 ± 11.2 90.5 ± 75.1 10.4 ± 7.8 -125.2 ± 73.5
SSP585 2030_2060 0 ± 0 9.8 ± 8.2 41.1 ± 18.4 19.1 ± 14.2 125.9 ± 134.4 12.2 ± 8.4 -164.8 ± 144.4
SSP245 2070_2100 0 ± 0 8.6 ± 11.2 44.2 ± 24.7 21.6 ± 19.9 123.1 ± 120.8 14 ± 9.3 -166.6 ± 139.5
SSP370 2070_2100 0 ± 0 19.2 ± 13.7 86.9 ± 44.6 37.7 ± 28 92.2 ± 110 30 ± 31.8 -177.7 ± 82.2
SSP585 2070_2100 0 ± 0 18.3 ± 19 81.8 ± 39.7 37.4 ± 24.7 150.7 ± 131 26.1 ± 24.6 -230.9 ± 148.4
$EUROPE-AFRICA
EUROPE-AFRICA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 -2.1 ± 5.2 -1.2 ± 9 -1.8 ± 19.9 1.8 ± 14.9 1 ± 20.7 -0.4 ± 8.2 13 ± 0
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 0.5 ± 5.8 0.2 ± 8.2 -1.4 ± 24.4 -4.2 ± 17.9 2.5 ± 16.7 2 ± 13.2 -6.5 ± 9.2
SSP245 2030_2060 11.2 ± 8.1 7 ± 10 26.5 ± 25.8 11.6 ± 19.8 -22.2 ± 35.5 -3.3 ± 10.1 -54 ± 0
SSP370 2030_2060 14 ± 6.8 6.1 ± 9.9 26.1 ± 20.3 10.5 ± 14.3 -23.1 ± 20.7 -4.5 ± 8 -56 ± 0
SSP585 2030_2060 16.9 ± 9.9 7.8 ± 10.2 35.5 ± 25.8 16.5 ± 17.2 -33.2 ± 32 -5.8 ± 13 -57 ± 0
SSP245 2070_2100 19.4 ± 8.2 6.8 ± 11.4 38.6 ± 29 14.9 ± 23.1 -36.5 ± 35.9 -2.5 ± 10.7 0 ± 0
SSP370 2070_2100 25.3 ± 16.7 -0.9 ± 23.1 35.5 ± 41.9 17.7 ± 22.1 -30.9 ± 28.5 -6.5 ± 21.3 0 ± 0
SSP585 2070_2100 32.4 ± 12.3 11.3 ± 19.9 56.3 ± 36.9 21.2 ± 22.1 -55.5 ± 27.4 -8.7 ± 18.1 0 ± 0
$NORTH-AMERICA
NORTH-AMERICA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 0 ± 0 0.4 ± 3.1 37 ± 71.9 31.8 ± 66.3 21.7 ± 91.3 4.9 ± 27.7 -56.1 ± 174
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 0 ± 0 0.3 ± 4.2 -5.4 ± 48.8 -3 ± 32.2 -19.8 ± 70 -2.7 ± 36.4 25.1 ± 114
SSP245 2030_2060 0 ± 0 3.4 ± 2.9 77.5 ± 109.8 60 ± 104.2 82.2 ± 94.5 14.2 ± 34.7 -159.2 ± 214.3
SSP370 2030_2060 0 ± 0 5.5 ± 6.9 62 ± 56.5 44.5 ± 45.2 70.5 ± 48.5 12.1 ± 33.2 -131.6 ± 77
SSP585 2030_2060 0 ± 0 7.1 ± 7.9 99 ± 102.3 62.1 ± 91 103.4 ± 97.9 29.8 ± 46.2 -201.3 ± 203.2
SSP245 2070_2100 0 ± 0 4.7 ± 3.1 87.1 ± 99.3 56.9 ± 86 145.9 ± 122.4 25.5 ± 49.7 -231.9 ± 215.6
SSP370 2070_2100 4 ± 0 15.7 ± 20 139.6 ± 80.2 91.5 ± 59.6 108 ± 110.3 28.4 ± 49.8 -243.9 ± 106
SSP585 2070_2100 5 ± 0 20.2 ± 38.4 147.1 ± 145.7 91.2 ± 136 198.8 ± 99.8 30.8 ± 35.4 -339.6 ± 220.2
$OCEANIA
OCEANIA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 1 ± 0 6 ± 42.8 4.5 ± 57.7 -1.2 ± 38.1 -3 ± 12.1 -1.2 ± 6.8 0.2 ± 1.2
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 1 ± 0 5.1 ± 32.9 5.6 ± 40.3 -2.2 ± 21.3 -3 ± 7.9 1.8 ± 9.4 -0.2 ± 1
SSP245 2030_2060 45 ± 0 50.2 ± 62.1 62.8 ± 74.9 -26.5 ± 39.9 -19.8 ± 18.3 -5.8 ± 12.4 -1.5 ± 0.7
SSP370 2030_2060 47 ± 0 42.4 ± 50.3 71.7 ± 61.2 -12.4 ± 33.3 -27 ± 61.2 -5.3 ± 10.3 -3.3 ± 3.2
SSP585 2030_2060 23.5 ± 0.7 33.8 ± 55.9 33.3 ± 71.6 -23.2 ± 43.1 -12 ± 20.1 -0.8 ± 12.2 -1.5 ± 0.7
SSP245 2070_2100 45 ± 0 66.6 ± 81.5 60.7 ± 101.2 -48.5 ± 58.5 -17.4 ± 18.8 -2.5 ± 13 -5 ± 0
SSP370 2070_2100 70.7 ± 84 63.2 ± 110.8 70 ± 153.2 -57 ± 62.9 -14.2 ± 57.5 -6.8 ± 13.5 -5 ± 0
SSP585 2070_2100 60.3 ± 48.7 69.8 ± 126.8 71.2 ± 164.5 -50 ± 91.8 -23.9 ± 22.1 -9 ± 14.2 0 ± 0
$SOUTH-AMERICA
SOUTH-AMERICA
model period Hyperarid Arid Sum drylands Semi-Arid Humid Dry subhumid Cold
historical 1850_1880 -0.1 ± 0.4 5.1 ± 13.9 26.3 ± 86.6 20.1 ± 83.4 -20.6 ± 79.5 1.3 ± 18.8 -6.8 ± 35.1
historical 1970_2000 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0 0 ± 0
historical 1985_2015 0.5 ± 1.2 2.7 ± 5.9 16.7 ± 19.2 10.5 ± 15.7 -14.2 ± 13.2 3 ± 9.3 -2.5 ± 5.8
SSP245 2030_2060 1.6 ± 1.1 13.2 ± 26.6 93.8 ± 100 55.6 ± 93.1 -78.5 ± 97.2 23.5 ± 24.8 -18.8 ± 42
SSP370 2030_2060 1.2 ± 0.9 5.9 ± 19.7 69.4 ± 52.4 38.2 ± 31.9 -59.5 ± 55.6 24.1 ± 36.6 -12.9 ± 18.7
SSP585 2030_2060 1.5 ± 1.1 19.2 ± 29.9 123.7 ± 97.4 72.8 ± 84.5 -107.5 ± 99.4 30.2 ± 38.1 -20.1 ± 41.2
SSP245 2070_2100 1.4 ± 1 22.6 ± 49.9 124 ± 111.5 68.9 ± 90.8 -107 ± 130.1 31 ± 41.3 -22.8 ± 43.7
SSP370 2070_2100 2 ± 2.3 22.9 ± 23.8 121.2 ± 92.7 73.5 ± 82.1 -109.5 ± 109.3 22.7 ± 35.9 -15 ± 15.1
SSP585 2070_2100 7.4 ± 15.3 40.5 ± 51.6 182.8 ± 140.1 96.9 ± 122.8 -160.2 ± 133.1 38 ± 40.7 -29.4 ± 47.8

10 Changes in proportion of aridity categories, figure

10.1 Mean changes of sum drylands by continent, period and model

list_plots <- list()

for(i in names(tab_list_percent)){
  df.percent <- tab_list_percent[[i]]
  df.sd <- tab_list_sd[[i]]
  
  df <- merge(df.percent, df.sd, by = c("model", "period"))
  
  g <- ggplot(data = df)+geom_point(aes(x=period,y= get("Sum drylands.x"), col = model), shape = "\u2605", size = 5, position = position_dodge(width = 0.5))+
    geom_errorbar(aes(x = period, ymin = get("Sum drylands.x")-get("Sum drylands.y"), ymax = get("Sum drylands.x")+get("Sum drylands.y"), col = model), position = position_dodge(width = 0.5))+
    scale_color_manual(values = colorvec[c(1,3,9,11)])+
    ylim(0,100)+
    labs(x="", y = "Percent of drylands", title = i)+
    theme_minimal()
  list_plots[[i]] <- g
}

ggarrange(plotlist = list_plots, ncol = 2, nrow = 4, common.legend = T, legend = "bottom")

10.2 Colplot by continent

tab.flow <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, Continent, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model, Continent) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL)

df245 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP245")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b245 <- ggplot(data = df245, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  facet_grid(rows = vars(Continent), switch = "y")+
  scale_y_continuous(position = "right")+
  labs(title = "SSP 2-4.5", y = "%", x = "")+
  theme_minimal()

df370 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP370")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b370 <- ggplot(data = df370, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  facet_grid(rows = vars(Continent), switch = "y")+
  scale_y_continuous(position = "right")+  labs(title = "SSP 3-7.0", y = "%", x = "")+
  theme_minimal()


df585 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP585")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b585 <- ggplot(data = df585, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  facet_grid(rows = vars(Continent), switch = "y")+
  scale_y_continuous(position = "right")+  labs(title = "SSP 5-8.5", y = "%", x = "")+
  theme_minimal()


ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)

10.3 Evolution of aridity index by Continent and period


tab.percent.region <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period != "1985_2015") %>%
  group_by(period, model, Continent, Name) %>% select(c("diff.AI.mean.percent")) %>%
  summarise(AI.diff.mean = mean(diff.AI.mean.percent, na.rm = T), AI.diff.sd= sd(diff.AI.mean.percent, na.rm = T)) %>%
  ungroup()

ipcc_regions.sf <- st_as_sf(ipcc_regions)

percent.region.sf <- merge(ipcc_regions.sf, tab.percent.region, all.y = T)
pbreaks <- c(-10, 0, 10, 20, 30, 40)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "azure","cornsilk", "#f4a582")
#colscale <- c("#c6dbef","white","#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#993344")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot(data = subset(percent.region.sf, period == j & model == i)) + geom_sf(aes(fill = AI.diff.mean))+
          binned_scale(aesthetics = "fill", breaks = pbreaks, palette = function(x) rev(colscale))+
      borders(colour = "grey60")+
      geom_sf_label(aes(label = round(AI.diff.mean,1)),size = 2) +
      labs(title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")


for(i in unique(percent.region.sf$Continent)){

  df <- percent.region.sf %>% filter(Continent == i & period != "1850_1880") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

  dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  
    
g <- ggplot(dfref)+geom_point(aes(x = period, y = AI.diff.mean, col = Name))+
  geom_line(aes(x = period, y = AI.diff.mean, col = Name, group = Name))+
  scale_color_manual(values = c(colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef","#e8ecff", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")))+
  facet_grid(cols = vars(model), scales = "free_x")+
  theme_minimal()+
  theme(legend.position = "right")
  
plot(g)
}

tab.percent.continent <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period != "1985_2015") %>%
  group_by(period, model, Continent) %>% select(c("diff.AI.mean.percent")) %>%
  summarise(AI.diff.mean = mean(diff.AI.mean.percent, na.rm = T), AI.diff.sd= sd(diff.AI.mean.percent, na.rm = T)) %>%
  ungroup()


df <- tab.percent.continent %>% filter(period != "1850_1880" & Continent != "POLAR") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  

dfref$Continent <- str_replace(dfref$Continent, "-", "\n")
    
ggplot(dfref, aes(x = period, y = AI.diff.mean))+
  geom_point(aes(size = AI.diff.sd, col = Continent))+
  geom_line(aes(group = Continent), size = 0.1)+
  facet_grid(rows = vars(model), cols = vars(Continent))+
  scale_size_continuous(breaks = c(0,1,5,10,100), transform = "pseudo_log", name = "Standard deviation")+
  scale_color_viridis_d(end = 0.8)+
  labs(y = "Difference in aridity index compared to 1970-2000, %", x = "")+
  guides(color = F)+
  ylim(-40,25)+
  theme_minimal()+
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 90))


tab.percent.continent <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period != "1985_2015") %>%
  group_by(period, model, Name) %>% select(c("diff.AI.mean.percent")) %>%
  summarise(AI.diff.mean = mean(diff.AI.mean.percent, na.rm = T), AI.diff.sd= sd(diff.AI.mean.percent, na.rm = T)) %>%
  ungroup()


df <- tab.percent.continent %>% filter(period != "1850_1880") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  

dfref$Name <- str_replace(dfref$Name, "-", "\n") 

ggplot(dfref, aes(x = period, y = AI.diff.mean))+
  geom_point(aes(size = AI.diff.sd, col = Name))+
  geom_line(aes(group = Name), size = 0.1)+
  facet_grid(cols = vars(model), rows = vars(Name), switch = "both")+
  scale_size_continuous(breaks = c(0,1,5,10,100), transform = "pseudo_log", name = "Standard deviation")+
  scale_color_viridis_d()+
labs(y = "Difference in aridity index compared to 1970-2000, %", x = "")+
  guides(color = F)+
  theme_minimal()+
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 90))+  
  scale_y_continuous(position = "right", limits = c(-40,30))

tab.percent.continent <- cmip6s %>% 
  subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period %in% c("1970_2000","2030_2060","2070_2100"))

df <- tab.percent.continent %>% filter(Continent != "POLAR") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  

dfref$Continent <- str_replace(dfref$Continent, "-", "\n")

ggplot(dfref, aes(x = period, y = diff.AI.mean.percent))+
  geom_boxplot(aes(col = Continent), outliers = F)+
  facet_grid(rows = vars(model), cols = vars(Continent), space = "free_y")+
  scale_color_viridis_d(end = 0.8)+
  labs(y = "Difference in aridity index compared to 1970-2000, %", x = "")+
  guides(color = F)+
  theme_minimal()+
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 90))

cmip6s$pr.percent <- with(cmip6s, diff.pr.mean/spr.ref.mean*100)

tab.percent.continent <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period != "1985_2015") %>%
  group_by(period, model, Continent) %>% select(c("diff.t.mean", "pr.percent")) %>%
  summarise(t.diff.mean = mean(diff.t.mean, na.rm = T), t.diff.sd= sd(diff.t.mean, na.rm = T),
            pr.diff.mean = mean(pr.percent, na.rm = T), pr.diff.sd= sd(pr.percent, na.rm = T)) %>%
  ungroup()


df <- tab.percent.continent %>% filter(period != "1850_1880" & Continent != "POLAR") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  

dfref$Continent <- str_replace(dfref$Continent, "-", "\n")
    
ggplot(dfref, aes(x = period))+
  geom_point(aes(y = t.diff.mean*3, size = t.diff.sd, col = "Temperature, °C"), alpha = 0.5)+
  geom_line(aes(y = t.diff.mean*3, group = Continent, col = "Temperature, °C"), size = 0.1)+
  
  geom_point(aes(y = pr.diff.mean, size = pr.diff.sd, col = "Precipitation, %"), alpha = 0.5)+
  geom_line(aes(y = pr.diff.mean, group = Continent, col = "Precipitation, %"), size = 0.1)+
  
  scale_color_manual(values = c("#3182bd", "#d6604d"))+
  facet_grid(rows = vars(model), cols = vars(Continent))+
  labs(y = "Precipitation anomalies, %", x = "", size = "Standard deviation", color = "")+
  ylim(-45,25)+
  theme_minimal()+
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 90))+
  
  scale_y_continuous(sec.axis = sec_axis(trans = ~./3, name = "Temperature anomalies, °C"))

11 Changes of variables over time and scenarios

11.1 Temperature

11.1.1 Absolute temperature

quantile(cmip6s$t.mean, probs = seq(0,1,0.1)) 
##          0%         10%         20%         30%         40%         50% 
## -53.4587313 -42.5828261 -28.4644864 -17.7502317  -7.2080554  -0.2580962 
##         60%         70%         80%         90%        100% 
##   7.1198406  15.3318920  22.8072968  26.4971129  34.3727997

t.breaks <- c(-50,-40,-30,-20,-10,0,10,20,30)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")
map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = t.mean))+
    borders(colour = "grey60")+
    binned_scale(aesthetics = "fill", breaks = t.breaks, palette = function(x) colscale,
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "Mean annual temperature, °C")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = t.mean))+
      borders(colour = "grey60")+
      binned_scale(aesthetics = "fill", breaks = t.breaks, palette = function(x) colscale,
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Mean annual temperature, °C", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}



ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")


#labels = c("2030-2060", "2070-2100"

11.1.2 Temperature anomalies

quantile(cmip6s$diff.t.mean, probs = seq(0,1,0.1)) 
##         0%        10%        20%        30%        40%        50%        60% 
## -3.5775507 -0.5616741  0.0000000  1.1670755  1.8544034  2.2566573  2.6600679 
##        70%        80%        90%       100% 
##  3.1444949  3.7454555  4.5875606 10.7627615

t.breaks <- c(-3,-0.5,0,1.5,2,3,4,5,10)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")
map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = diff.t.mean))+
    borders(colour = "grey60")+
   scale_fill_gradient2(low = "#08519c",mid = "white", high = "#b2182b", midpoint = 0, transform = "pseudo_log", breaks = c(-2, -1, 0, 1, 2, 3, 4))+
      labs(fill = "Precipitation anomalies in %", title = i)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = diff.t.mean))+
      borders(colour = "grey60")+
   scale_fill_gradient2(low = "#08519c",mid = "white", high = "#b2182b", midpoint = 0, transform = "pseudo_log", breaks = c(-2, -1, 0, 1, 2, 3, 4))+
      labs(fill = "Temperature anomalies in °C", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

### Binned anomalies

YOR <- brewer.pal(n = 11, name ="YlOrRd")
coltemp <- c("[0,1]" = YOR[1], "(1,2]" = YOR[2], "(2,3]" = YOR[3], "(3,4]" = YOR[4], "(4,5]" = YOR[5], "(5,6]" = YOR[6], "(6,7]" = YOR[7], "(7,8]" = YOR[8], "(8,9]" = YOR[9], "(9,10]" = "grey20")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot(data = subset(cmip6s, period == j & model == i), 
                aes(x=lon, y = lat,  fill = cut(diff.t.mean, breaks = c(0,1,2,3,4,5,6,7,8,9,10), include.lowest = T))) + 
      geom_tile(show.legend = T)+
      borders(colour = "grey60")+
      scale_fill_manual(values = coltemp, drop = F)+
      labs(fill = "Temperature\nanomalies in °C", title = index)+
      theme_void()+ylim(-55,90)+
     theme(legend.position = "bottom", legend.text =  element_text(angle = 0))+
    guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

11.2 Precipitation

11.2.1 Absolute values of precipitations

quantile(cmip6s$pr.mean, probs = seq(0,1,0.1)) 
##          0%         10%         20%         30%         40%         50% 
##    8.579696   55.000319  104.132310  200.852004  337.042710  463.946055 
##         60%         70%         80%         90%        100% 
##  579.741473  712.099280  959.490392 1431.315602 6863.723388

p.breaks <- c(10,50,100,200,300,400,500,600,700)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")
map_list <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = pr.mean))+
    borders(colour = "grey60")+
    binned_scale(aesthetics = "fill", breaks = p.breaks, palette = function(x) rev(colscale),
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "Annual mean\nprecipitation, mm")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = pr.mean))+
      borders(colour = "grey60")+
      binned_scale(aesthetics = "fill", breaks = p.breaks, palette = function(x) rev(colscale),
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Annual mean\nprecipitation, mm", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

11.2.2 Precipitation anomalies

cmip6s$pr.percent <- with(cmip6s, diff.pr.mean/spr.ref.mean*100)

quantile(cmip6s$pr.percent, probs = seq(0,1,0.1)) 
##           0%          10%          20%          30%          40%          50% 
## -42.02999493  -2.87352237   0.00000000   0.04663897   2.46136429   5.19906939 
##          60%          70%          80%          90%         100% 
##   9.34748067  13.67239987  18.11758968  24.02085190 148.15226168

p.breaks <- c(-43,-3,0,3,5,9,13,18,24,148)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")
map_list <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = pr.percent))+
    borders(colour = "grey60")+
   scale_fill_gradient2(high = "#08519c", mid = "white", low = "#b2182b", midpoint = 0, transform = "pseudo_log", breaks = c(-10, -5, -1, 0, 1, 5, 10, 15, 30))+
      labs(fill = "Precipitation anomalies in %", title = i)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = pr.percent))+
      borders(colour = "grey60")+
   scale_fill_gradient2(high = "#08519c", mid = "white", low = "#b2182b", midpoint = 0, transform = "pseudo_log", breaks = c(-10, -5, -1, 0, 1, 5, 10, 15, 30))+
      labs(fill = "Precipitation anomalies in %", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = pr.percent))+
      borders(colour = "grey60")+
   scale_fill_binned_divergingx(mid = 0, palette = "RdBu", rev = F, n_interp = 7, breaks = c(-20, -10, 0, 10, 20, 30, 40, 50, 60))+
      labs(fill = "Precipitation anomalies in %", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")